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ABSTRACT 


A procedure is presented to perform a contact analysis of spiral bevel gears in order to 
predict the contact path and the load sharing as the gears roll through mesh. The approach 
utilizes recent advances in automated contact methods for nonlinear finite element analysis. 
A sector of the pinion and gear is modeled consisting of three pinion teeth and four gear 
teeth in mesh. Calculation of the contact force and stresses through the gear meshing cycle 
are demonstrated. Summary of the results are presented using 3-dimensional plots and 
tables. Issues relating to solution convergence and requirements for running large finite 
element analysis on a supercomputer are discussed. 
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CHAPTER I 


INTRODUCTION 

Spiral bevel gears are important elements for transmitting power between intersecting shafts. 
They are commonly used in applications that require high load capacity at high operating speeds. 
One such application is in helicopter transmission systems. Aircraft designers are continually 
required to improve performance. Reduced weight, size, and cost with increased power, 
efficiency, and reliability are constantly being sought. 

Prior research has focussed on various aspects of spiral bevel gear operation. Much has 
been done on spiral bevel gear geometry to reduce noise and vibration, kinematic error, improve 
manufacturability, and inspection. Stress analysis is another important area of ongoing 
research. Accurate prediction of contact stresses and tooth root/fillet stresses are important to 
increase reliability and reduce weight. 

Much effort has focussed on predicting stresses in gears with the finite element method. 
Most of this work has involved parallel axis gears with two dimensional models. Only a few 
researchers have investigated finite element analysis of spiral bevel gears (ref. 1-4). In reference 
4, finite element analysis was done on a single spiral bevel gear tooth using an assumed contact 
stress distribution. In that model, contact stresses were not evaluated. 

For parallel axis gears, a closed form solution exists which determines the surface 
coordinates of a tooth. This is then used as input to finite element methods. For spiral bevel 
gears there is no closed form solution. Therefore, the kinematics of the cutting or grinding 
machinery is utilized to numerically describe the surface coordinates of the gear tooth. 

The research reported herein is based on the extension of work presented in references 4-7. 
A model that has three pinion teeth and four gear teeth has been developed based on gear 
manufacturing kinematics for a single tooth on each the pinion and the gear. 
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The objective of this research is to study the contact path and load sharing in these gears 
when contact occurs on multiple teeth in mesh. This is done by performing a static analysis at 
different incremental rotations. A nonlinear approach is required due to large displacements 
associated with gear rotation and nonlinear boundary conditions associated with the gear tooth 
surface contact. Also evaluated are the contact stresses and fillet stresses. 
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CHAPTER II 


GEAR SURFACE GEOMETRY 

Briefly described in this chapter are the gear manufacturing process, the kinematics of the 
manufacturing process, tooth surface coordinates solution technique, surface rotations of the gear 
and pinion, and different orientations required for the spiral bevel gears to mesh with each other. 

2.1 Gear Manufacture 

Machinery for the manufacture of spiral bevel gears is provided by the Gleason Works, 
Rochester, NY. These machines are preferred over gear hobbing machines because they can 
be used for both milling and grinding operations. Grinding is important in producing hardened 
high quality aircraft gears. 

This machine consists mainly of three parts: the machine frame, the cradle, and the sliding 
base (ref. 8). The cradle with the head cutter mounted on it, slowly rotates about its axis, as 
does the gear which is being cut. During this motion the gear surface is generated. As the 
cutter disengages from the workpiece, the cradle reverses rapidly and the sliding base translates 
with respect to the cutter to index the workpiece for the next cutting cycle. This sequence is 
repeated until the last tooth is cut. 

The head cutter is mounted on the cradle with an offset from the cradle center. This allows 
adjustment of the axial distance between the cutter center and the machine center. The 
adjustment of the angular position between the two axes provides the desired spiral angle. The 
shape of the blades of the head cutter are typically straight lines that rotate about their own axis 
at a speed for efficient metal removal. The rotational speed of the cutting head is independent 
of the cradle or workpiece rotations (ref. 9). 
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The pinion is typically cut one side at a time, whereas the gear is cut both sides 
simultaneously. Spiral bevel gears can be either left or right handed. In a left handed spiral 
bevel gear, the tooth spirals to the left while looking from the apex of cone towards the gear. 
Whereas, for a right handed spiral bevel gear the tooth spirals to the right. 


2.2 Tooth Surface Coordinate Equations 

The system of equations, required to define the coordinates of spiral bevel gear tooth 
surfaces, were derived in reference 4, and are briefly summarized here. 

The first equation is the equation of meshing. This equation is based on the kinematics of 
manufacture and the machine tool settings. The equation of meshing requires that the relative 
velocity between a point on the cutting surface and the same point on the pinion being cut must 
be perpendicular to the cutting surface normal (ref. 9). 

n • V = 0 

where, n is the normal vector from the cutter surface and 

V is the relative velocity between the cutter and the workpiece surfaces 
at the specified location. 

This equation is developed in terms of machine tool coordinates U, © and <f> e 

r cotilf - U cos\|r 
U sini|r sin0 
c U sinijr cos0 
1 

where, U is the generating cone surface coordinate used to locate a point along the length 

of the cutting head 

0 is the generating cone surface coordinate used for rotational orientation of a 
point on the cutting head 

<f> t is the rotated orientation of the cutter as it swings on the cradle 
r is the radius of the generating cone surface and is described by the following 
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equations (ref. 4): 

Since the kinematic motion of cutting a gear is equivalent to the cutting head meshing with a 
simulated crown gear, an equation of meshing can be written in terms of a point on the cutting 
head i.e in terms of U, 0 and </> e . The equation of meshing for straight sided cutters with a 
constant ratio of roll between the cutter and the workpiece is given by (ref. 1): 


(U - i coti|f cosi|r) cosy sinx 
+ - siny) cosi|r sin0 

* cosy siny sin(g- - 4> c ) (2) 

± E m (cosy siny + siny cosy cost) ' 

- L m siny cosij; sint = 0 

The upper and lower signs are for left and right hand gears respectively. The following 
machine tool settings are defined: 


r (O-Fqi^c) 

q cradle angle 

7 root angle of workpiece 

Ejj machining offset 

L m vector sum of change of machine center to back and the sliding base 
m cw = <£ c /<£ w , the relationship between the cradle and the workpiece for a constant ratio 
of roll 

U generating cone surface coordinate 

S radial location of cutting head in coordinate system S m 

r radius of generating cone surface 


Equation (2) is equivalent to: 


f x {U, 6, <{> c ) = 0 


(3) 


Because there are three unknowns U, 0, and </> c ; three equations must be developed to solve 
for the surface coordinates of a spiral bevel gear. The three parameters U, 0 and <f> c are 
defined relative to the cutting head and cradle coordinate systems (S c and S s ) respectively. 
These parameters can be transformed through a series of coordinate transformations to a 
coordinate system attached to the workpiece. Or U, 0 and 0 C can be mapped into X w , Y w , Z w 
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in coordinate system S w attached to the workpiece. These transformations, used in conjunction 
with two other geometric requirements, give the two additional equations. 

The correct U, 0 and <f> c that solves the equation of meshing, must also, upon transformation 
to the workpiece coordinate system S w , result in a axial coordinate Z w that matches with the 
value of Z found by the projection of the tooth in the XY plane. 

Z w - Z = 0 ( 4 ) 

This equation along with the correct coordinate transformations (see Equation 1 1) result in a 
second equation of the form: 

f 2 (u, e, <j> c ) =o (5) 

A similar requirement for the radial location of a point on the workpiece results in the 
following: 

I - + Yl) = o ( g ) 

This is shown in figure 2.1. The appropriate coordinate transformations (see Equation 1 1) 
will convert equation 6 into a function of U, 0 and 4> c 

f 3 (u, e, 4> c ) = o (7) 

Equations (3), (5), (7) form a system of nonlinear equations necessary to define a point on 
the tooth surface. 
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2.3 Solution Technique 

The three equations discussed earlier to describe the tooth surface coordinates are nonlinear 
and do not have a closed form solution. They are solved using Newton’s method (ref. 10). 

An initial guess U 0 , 0 O , 4> c0 is used to the start iterative solution procedure. Newton’s 
method is used to determine subsequent values of the updated vector (U k , 0 k , </> ck ). 


U k 




Yk-i 


= 


+ 

Yk-i 





xU 


where the vector Y is the solution of: 


( 8 ) 


dfiUr*- 1 ) 


0f 2 (0*- 1 ) 

dU 

06 



df 2 (6*" 1 ) 


dU 

00 


df 3 ( U k_1 ) 

dfj (0* -1 ) 

df 3 (Of 1 ) 

dU 

00 

0<J) C 





if' 1 


f x ( U*~ l , 6* -1 , 4>e 1 ) 

r 2 

► = < 

£■% ( u k_1 , 0* _1 , 4>c _1 ) < 

if" 1 




The 3x3 matrix in the preceding equation is the Jacobian matrix and must be inverted each 
iteration to solve for the Y vector. The equation of meshing, function fj, is numerically 
differentiated directly to find the terms for the Jacobian matrix. Function f 2 and f 3 cannot be 
directly differentiated with respect to U, 0 and <£ c . After each iteration U*' 1 , 0 k ‘‘, (in the 
cutting head coordinate system) are transformed into the workpiece coordinate system) and are 
transformed into the workpiece coordinate system, S w , with the series of coordinate 
transformations as given in Equation 1 1 . 
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( 10 ) 


x w 



Y 


r coti|r - U cosi|r 

1 w 
rr 

= t«wc] 

U sint|r sin6 



U sini|r cos0 

.1. 


** 


where, 

[M wc ] « [M wa f (4> e ) [M ap 3 [j^J (4) c ) ] [M-SC] (11) 

Each matrix [M] above represents a transformation from one coordinate system to another 
(ref. 4). 

Functions f 2 and f 3 are evaluated by starting with an initial U k , 0 k and <f> k , performing the 
transformations in Equation (11) and evaluating Equations (4) and (6). The numerical 
differentiation of f 2 and f 3 is performed by transforming U k +inc, 0 k +inc, 0 k +inc (where inc 
is a small increment appropriate for numerical differentiation) into ^4- inc, Y w +inc, Z w +inc. 
Equations (4) and (6) are then used to evaluate the numerical differentiation. Function fj, f 2 , 
f 3 and their partial derivatives are required for the Jacobian matrix and are updated each 
iteration. The iteration procedure continues until the Y vector is less then a predetermined 
tolerance. This completes the solution technique for a single point on the spiral bevel gear 
tooth surface. In this way the coordinates of the surface of the tooth are defined. This 
solution technique is repeated four times for each of the four surfaces; gear convex, gear 
concave, pinion convex and pinion concave. 

Since all four surfaces are generated independently, additional matrix transformations are 
required to obtain the correct tooth thickness. The concave surfaces are fixed on each tooth. 
The convex surfaces are rotated as required. The angle of rotation is obtained by matching the 
tooth top land thickness with the desired value. 
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2.4 Gear and Pinion Orientations Required for Meshing 

After generating the pinion and gear surface as described above, the pinion cone and gear 
cone apex will meet at the same point (see figure 2.2). This point is the origin of the fixed 
coordinate system attached to the workpiece being generated. To place the gear and pinion in 
mesh with each other, rotations described in the following example are required (ref. 5). 

1. The gear tooth surface points are rotated by 360/N t +180 CW about the global 7 ^ axis, 
for this example, the rotation is 190 degrees. 

2. The pinion is rotated by 90 deg CCW about the global Y axis. 

3. Because the four surfaces are defined independently, their orientation is random with 
respect to meshing. The physical location of the gear and pinion after rotations described 
above correspond to the gear and pinion overlapping. To correct this condition the 
pinion is rotated CW about its axis of rotation Z w until surface contact occurs. For the 
example used in this study, the rotation was 3.56 deg. 

To make a complete pinion, the generated pinion tooth was copied and rotated 12 times and 
the generated gear tooth was copied and rotated 36 times with the aid of a solid modelling 
program as shown in figure 2.3. 
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CHAPTER III 


CONTACT ANALYSIS BY THE FINITE ELEMENT METHOD 

The advantages of the Finite Element Analysis (FEA) for accurate deformation and stress 
analysis of complex forms is well known. This chapter provides a brief outline of the FEA 
analysis carried out for spiral bevel gears. It also describes the fundamental concepts of non- 
linearity with emphasis on automated contact analysis. Since the research reported herein is 
presented using the general purpose finite element code MARC (ref. 11), details of its special 
features used for the analysis and the description of various blocks of the input deck are also 
discussed. 

3.1 FEA of Spiral Bevel Gears 

Only a few researchers have investigated finite element analysis of spiral bevel gears. FEA 
analysis has been done on a single spiral bevel tooth using an assumed contact stress distribution 
(ref. 4). More recent FEA spiral bevel gear analysis research has attempted to solve the contact 
stress distribution in a multi-tooth model (4 gear and 3 pinion teeth) (ref. 3). The tooth pair 
contact zones in reference 3 were modeled with gap elements. The study here uses software 
with automated contact options in order to avoid certain limitations in the use of gap elements, 
such as: 

(i) It is difficult to connect the gap elements with proper orientations across the 
normal from one surface to the other surface parallel to the contact point. 

(ii) The orientation of the contact plane remains unchanged during deflection. 

(iii) It is difficult to accurately select the properties such as appropriate open/closed 
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stiffness values, selection of the stiffness matrix update strategy and efficient 
problem restarts. 

New advances in the state of the art for FEA provide deformable body against deformable 
body penetration algorithms which can be used to establish the nonlinear boundary conditions 
for contact problems. MARC (ref. 1 1) is one such FEA package software which uses this 
algorithm to automatically detect contact. 

3.2 Concepts of Nonlinear Analysis 

In linear FEA, a simple linear relationship exists between force and deflection (Hooke’s 
law). For a metallic spring under small strain, the force F equals the product of the stiffness 
K and the deflection U. 

Also, the deflection can be obtained by dividing the force with the spring stiffness. This 
relationship is valid as long as the spring remains linear elastic, and the deflections are such that 
they do not cause the spring to yield or break. If the spring material is changed, for example, 
from steel to rubber, the linear force- displacement relation is no longer valid. It becomes a 
nonlinear problem i.e. F=£KU. 

A nonlinear structure is the one for which the force-deflection relationship cannot be directly 
expressed in terms of a set of linear equations. 

The three major types of nonlinearities are: 

(i) Geometric nonlinearity (large deformations, large strains, snap through buckling) 

(ii) Material nonlinearity (plasticity, creep, viscoelasticity) 

(iii) Boundary nonlinearity i.e., a changing status (opening/closing of gaps, contact, 
follower force) 

A nonlinear system cannot be analyzed directly with a linear equation solver. However, 
it can be analyzed by using a series of linear approximations. Each linear approximation 
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requires a separate pass or iteration, through the program’s linear equation solver. Each new 
iteration is as expensive, in terms of CPU time, as a single linear analysis solution. 

In the preprocessing phase of a nonlinear analysis, everything is quite similar to linear FEA 
data input except the user is required to specify certain nonlinear analysis controls (i.e., large 
displacements, "contact" parameters, convergence controls, etc.) and additional material 
properties required for a nonlinear analysis. 

In the solution phase, of nonlinear FEA, the solver must perform the analysis in steps or 
increments. Within each increment, the code will also iterate as necessary until equilibrium is 
achieved, before proceeding to the next increment. 

In the post processing phase, the user looks at quantities like stress contours, etc. A 
nonlinear analysis takes tens, hundreds and sometimes, thousands of increments, thus, usually 
requiring a high speed computer with lots of memory for a reasonable turnaround. The 
objective in a successful nonlinear analysis is to obtain a converged solution at a reasonable cost. 

In large deformation analysis, incremental load AF and displacement AU is related by the 
tangent stiffness K x . In solving this type of problem, the load is increased in small increments, 
the incremental displacement AU is found, and the next value of the tangent stiffness is found 
and so on. A brief description of the incremental solutions will now be discussed. 

FEA is an approximate technique, and there exist many methods to solve the basic equations. 
In nonlinear FEA, two popular incremental methods used to solve the nonlinear equilibrium 
equations are: Full Newton-Raphson or the Modified Newton-Raphson. 

The Full Newton-Raphson Method (see figure 3.1) assembles and solves the stiffness matrix 
at every iteration, and is thus expensive for large 3-D problems. It has quadratic convergence 
properties, which means in subsequent iterations the relative error decreases quadratically. It 
gives good results for most nonlinear problems. 

The Newton-Raphson method iterates using this equation: 
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[K t ] {AC/} = {F^} - {FJ 


(ID 


where, 

[K t ] = the tangent stiffness matrix 

{AU} = the displacement increment 

{F app } = the applied nodal force 

{For} = the restoring N-R force (the loads generated by the elemental stresses) 

{F app }-{F 0 } = the residual force 

The program updates the tangent stiffness matrix [K x ] and the residual F app -F nr at each 
iteration, and then resolves the equation given above. Convergence is achieved once F app -F nr is 
less then a convergence criterion that is set. If F app is not equal to F^, the system is not in 
equilibrium. 

As shown in figure 3.1, the 1st iteration yields a displacement AU„ using the initial elastic 
stiffness and the applied load F app . The nonlinear response yields a force value F„ for this 
displacement. The 2nd iteration yields AU 2 , using the updated tangent matrix and the residual 
load. Subsequent iterations quickly drive the analysis to a convergent solution. This solution 
guarantees convergence if, and only if, the starting AU is near the exact solution. This could 
be obtained by taking smaller load increments. 

In solving the convergence problem, one cannot neglect the general FEA conflict of expense 
versus accuracy. One must balance computational expense against accuracy. Using a finer 
mesh and multiple load increments can often lead to a more accurate and more robust (less likely 
to diverge) solution, but usually at increased expense. 

The Modified Newton-Raphson method does not reassemble the stiffness matrix during every 
iteration as shown in figure 3.2. It costs less per iteration, but the number of iterations may 
increase substantially over that of the full Newton-Raphson method. It is effective for mildly 
nonlinear problems. In this analysis the full Newton-Raphson method was used. 
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3.3 Automated Contact Analysis 

Many common structural features exhibit nonlinear behavior that is status-dependent. The 
stiffness of these systems shifts abruptly between different values, depending on the overall 
status of the item. Status changes might be directly related to load, or they might be 
determined by some external cause. 

Situations in which contact occurs are common to many different nonlinear applications. 
Contact forms a distinctive and important subset to the category of changing-status nonlinearities 
(ref. 15). Contact, by its very nature, is a nonlinear problem. During contact, both the forces 
transmitted across the area and the area of contact change. The force-displacement relationship 
thus become nonlinear. Usually, the transmitted load is in the normal direction. In this report 
the method of reference 1 1 is used to perform the nonlinear analysis. 

Reference 1 1 is a general purpose computer program for linear and nonlinear stress and heat 
transfer analysis. This program is capable of solving problems with nonlinearities that occur due 
to material properties, large deformations, or boundary conditions. In general, the solution of 
nonlinear FEA problems requires incremental solution schemes and sometimes requires iterations 
within each load/time increment to satisfy equilibrium conditions at the end of each step. The 
program features relevant to gear meshing are discussed in this section. 

The FEA program used has a fully automatic CONTACT option which enables the analysis 
between finite element bodies without the use of any special gap or contact elements. The 
procedure is capable of handling the following types of contact: 

i) deformable bodies against rigid surfaces 

ii) deformable bodies against deformable bodies 

iii) a deformable body against itself 

The CONTACT option was originally designed for analysis of manufacturing processes such as 
forging or sheet metal forming, but its capabilities have been expanded to meet other analysis 
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requirements. The work presented in this document utilized the program of reference 11 
running on a supercomputer. 

Contact between the bodies is handled by imposing non-penetration constraints (reference 
11). The non-penetration constraint, as shown in figure 3.3 and figure 3.4, is 

U A • n < D 

A non-penetration constraint (that the surfaces cannot inter-penetrate) is usually handled in FEA 
codes by one of three techniques: Lagrange multipliers (used in several FEA codes that offer 
"gap" elements); penalty functions (one application being the use of stiff or rigid connecting 
members to approximate the constraint); or solver constraints. In some FEA codes which offer 
explicit dynamics capabilities, a fourth technique is the direct application of contact forces. Use 
of a "gap element" means node to node contact. The CONTACT features uses the solver 
constraints approach. 

Solver constraints are used to impose the non-penetration constraint, and a very efficient 
surface contact algorithm which allows the user to simulate general 2 and 3D multibody contact. 
Both "deformable-to-rigid" and deformable-to-deformable" contact situations are allowed. The 
user needs to only identify which bodies might come in contact during the analysis. Self- 
contact, which is common in rubber problems, is permitted. The bodies can be either rigid or 
deformable, and the algorithm tracks variable contact conditions automatically. Thus the user 
no longer needs to worry about the location and open/close status checks of "gap elements". 
Automatic, in this context, means that user interaction is not required in treating multibody 
contact and friction, and the program has automated the imposition of non-penetration 
constraints. Also, coupled thermo-mechanical contact problems (e.g rolling, casting, forging) 
and dynamic contact problems can be handled. 

Real world contact problems between rigid and/or deformable bodies are actually 3D in 
nature. To solve such contact problems, one needs to define bodies and their boundary 
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surfaces. The definition of bodies is the key concept in analyzing 3D contact automatically. 
For rigid bodies, one can define such surfaces as: 4-point patch; ruled surface; plane; tabulated 
cylinder; surface of revolution; and Bezier surfaces. 

Deformable bodies are defined by the elements of which they are made. Once all the 
boundary nodes for a deformable body are determined, 4 point patches are automatically created, 
which are constantly updated with the body deformation. Contact is determined between a node 
and all body profiles, deformable or rigid. A body may fold upon itself, but the contact will 
still be automatically detected, thus preventing self penetration. 

The user must define bodies so that their boundary surfaces can be established. Deformable 
bodies are defined by a list of finite elements, and rigid bodies are defined by geometrical 
entities. Because the contact boundary conditions are a function of the applied load, the 
analysis must be carried out incrementally. Within each load or time increment of an analysis, 
additional iterations may be required to stabilize the contact zone. Contact problems involve 
two important aspects: 

i) the opening and closing of the gap between bodies 

ii) friction between the contacting surfaces. 

The MARC program establishes a hierarchy between the bodies so that at a given contact 
interface, one body is the contactor and the other body is the contacted. The set of nodes on 
the boundary surface of a contactor are candidate nodes for contact. The boundary surface of 
a contacted body is defined by a set of geometrical entities. A user specified contact tolerance 
is used to determine the body separation distance which determines whether two bodies are 
considered to be in contact with each other. The contactor’s boundary nodes are prevented 
from penetrating the surface of the contacted body by imposing solver constraints. For contact 
between a deformable body and a rigid body, a displacement constraint is applied. For contact 
between two deformable bodies MARC applies multipoint constraints in the form of ties. The 
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ties link the motion of one node on the contactor body to two adjacent nodes on the contacted 
body. During each iteration as nodes enter and leave contact or slide between adjacent node 
segments on contacted bodies, a bandwidth optimization is performed to reduce the computer 
processing time required. 

During contact, it is possible to move bodies around during an analysis; however, the user 
must make sure that deformable bodies do not have rigid body modes. A minimum number 
of boundary conditions or spring element attachments must be applied to prevent rigid body 
motion. 

A static analysis of two (or) more bodies that are not initially connected poses special 
problems with a finite element analysis, if one of the bodies has a rigid body motion component. 
If, at any time, the two bodies are disconnected then the stiffness matrix would become singular 
and unsolvable (in a static analyses). This is because finite elements require at least some 
stiffness connecting all the elements together along with sufficient displacement constraints to 
prevent rigid body motion. In order to overcome this difficulty, weak springs have been added 
to connect the bodies. The spring stiffness is very low. This stiffness is there only to supply 
some stiffness to the system. The stiffness should be negligible compared to the material 
stiffness, so that it has no effect on the solution. 
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CHAPTER IV 


FINITE ELEMENT MODEL AND ANALYSIS 

This chapter describes the procedure for assembling a finite element gear pair model for 
analysis with reference 1 1 software. Two different models are analyzed. These models are 
generated in the geometric modelling program of reference 14 (FEA pre and post-processor). 
The first model is a two tooth segment and the second model is a seven tooth segment of a 
gearset. Different programs used to create these FEA models and the boundary conditions 
imposed on them are described. Various sets of rotations carried out in the gearsets for the 
analysis to determine the contact path and the load sharing across the tooth are then discussed. 
A detailed report of the problems faced during the course of the analysis for convergence to take 
place and memory and CPU requirements of the system used are described. Different 
parameters of the finite element analysis input commands, which have been studied to reduce 
the CPU and memory requirements of the system, are also reported in this chapter. 

4.1 Model Descriptions 

To model spiral bevel gears, the machine settings and the gear and pinion design data as 
given in Tables I and II are necessary. The equation of meshing and the kinematics of gear 
cutting are incorporated in the computer programs to generate the spiral bevel gear model in the 
geometric modeling program (ref 14). The input data (for the geometric modeling program) 
for the seven tooth model is obtained from references 16, 17. This is a 10 x 10 mesh input to 
generate 4 gear teeth and 3 pinion teeth in mesh to simulate contact on multiple teeth of a spiral 
bevel gearset. The input also includes the hub attached to the gear. The input data (for the 
geometric modeling program) for the two tooth model is obtained from reference 6 and 7. The 


18 



programs generate a NxM mesh on the tooth profile of a spiral bevel gearset. 


TABLE I: PINION AND GEAR DESIGN DATA 


PINION 

GEAR 

Number of teeth pinion 

12 

36 

Dedundum angle, degree 

1.5666 

3.8833 

Addendum angle, degree 

3.8833 

1.5666 

Pitch angle, degree 

18.4333 

71.566 

Shaft angle, degree 

90.0 

90.0 

Mean spiral angle, degree 

35.0 

35.0 

Face width, mm (in) 

25.4(1.0) 

25.4 (1.0) 

Mean cone distance, mm (in) 

81.05 (3.19) 

81.05 (3.19) 

Inside radius of blank, mm (in) 

5.3 (0.6094) 

3.0 (.3449) 

Top land thickness, mm (in) 

2.032 (0.080) 

2.489 (0.098) 

Clearance, mm (in) 

0.762 (0.030) 

0.92964 (0.0366) 


TABLE II: GENERATION MACHINE SETTINGS 



CONCAVE 

CONVEX 

CONCAVE 

CONVEX 

Radius of cutter, r, in 

2.9656 

3.0713 

3.0325 

2.9675 

Blade angle, i/s degree 

161.9543 

24.33741 

58.0 

22.0 

Vector sum, L«, 

0.038499 

-0.051814 

0.0 

0.0 

Machine offset, E m 

0.154576 

-0.1742616 

0.0 

0.0 

Cradle to cutter distance, s, in 

2.947802 

2.8010495 

2.285995 

2.285995 

Cradle angle, q, degree 

63.94 

53.926 

59.234203 

59.234203 

Ratio of roll, M^ 

0.30838513 

0.32204285 

0.950864 

0.950864 

Initial cutter length, u, in 

9.59703 

7.42534 

8.12602 

7.89156 

Initial cutter orientation, 6, degree 

126.83544 

124.43689 

223.9899 

234.9545 

Initial cutter orientation, 4> c , degree 

-0.85813 

-11.38663 

-0.35063 

-12.3384 
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A 24 x 12 mesh was used to create the two tooth model. Eight noded, isoparametric HEX 
elements were used. The seven tooth model and the two tooth model are as shown in figures 

4.1 and 4.2. The seven tooth model consisted of 8793 elements and 11261 nodes and the two 
tooth model consisted of 3116 elements and 4452 nodes. 

4.2 Loading and Boundary Conditions 

The boundary conditions for the seven teeth and two teeth models are as shown in figures 
4.1 and 4.2. The boundary conditions are applied such that the gear teeth are made to pivot 
about a fixed point at node number 7872 for the seven tooth model and at a node number 4448 
for the two tooth model. In both models, the axis of rotation for the gear is the Z axis. The 
nodes where the forces are applied are in the gear hub of the models. Fixed displacement 
boundary conditions are applied at 8 nodal points in the Z direction only for the gear and in all 
directions for the pinion as shown in the figures. Since this is a static problem involving two 
bodies (the pinion and the gear) in contact, as described in a previous chapter, weak springs are 
added to prevent the rotational rigid body modes for the gear. Eight springs are added away 
from the contact region. The springs connect the comer nodes of the pinion with the comer 
nodes of the gear on the faces where contact occurs. The stiffness of the springs are 100 lbs/in. 

This is insignificant when compared to the tooth contact stiffness and therefore does not effect 
the overall solution. 

The maximum torque for the gear mesh studied was 9508 in.» lbs. on the gear. This torque 
is applied as a concentrated force with a moment arm on the gear hub. This concentrated force 
for the seven tooth model was 4714 lbs. with a moment arm of 2.017 inches. For the two tooth 
model, the force is 3392 lbs. applied at a radius of 2.798 inches. The force is applied in 
different increments for convergence to occur. The details of the incremental loading are 
discussed later in this chapter. 
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4.3 Model Generation to Predict the Contact Path 

To predict the contact path of the spiral bevel gears as they roll in and out of mesh, several 
rotations have been carried out on the preliminary model. The rotation of the gear and the 
pinion should be in the same directions, i.e., both positive. 1 For the model being analyzed, 
the gear has 36 teeth and the pinion has 12 teeth. Therefore, for each degree of gear rotation 
the pinion has to be rotate three degrees to avoid interference during meshing. For the seven 
tooth model, seven such rotations were carried out. For the two tooth model, three such 
rotations have been carried out. From the positioning of the preliminary model, the rotations 
were carried out in both positive and negative directions to determine the contact path. In the 
seven tooth model the pinion was rotated by +6, 0,-6, -12, -18, -24, -30 degrees and in the two 
tooth model the pinion was rotated by +6, 0,-6 degrees. The corresponding rotation of the gear 
were -2,0, +2, +4, +6, +8, + 10 degrees for the seven tooth model and -2,0, +2 degrees for the 
two tooth model. (Zero degrees corresponds to the initial position for the gearset based on 
solving the equation of meshing.) 

While generating the models, care must be taken not to change the node numbers and the 
element numbers. This eliminates a great deal of editing in the FEA input file. The command 
used in the geometric modeling program to do this with the gear and pinion part names is as 
follows: 

Name, (partname), rot, rotation axis data, (partname) 

Care must also be taken to rotate the gear about the Z axis and pinion about the X axis 
directions to avoid interference of the gears in mesh. 


'NOTE: Although all gearsets must rotate in the opposite sense, (i.e., if the gear rotates CW, the pinion must 
rotate CCW) this gearset does in fact rotate in the same direction (i.e., both gear and pinion have positive rotation) 
based on the right hand rule and the defined coordinate system. This occurs because the gear rotates about the Z 
axis, with the positive Z axis pointing from the toe of the gear to the heel. However, the pinion rotates about the 
X axis with the positive X axis pointing from the heel of the pinion to the toe. The result is both gear and pinion 
have positive rotation. 



4.4 Assembling the Spiral Bevel Gear Pair for FEA Analysis 

The procedure to assemble a spiral bevel gear pair to perform FEA analysis is as follows: 
STEP 1. Initially the gear model is generated in the geometric modelling program, from the 
input file obtained after executing the computer codes from references 6 and 7 as discussed 
earlier. The model is optimized and the node and the element id’s are compacted. The load 
and the boundary conditions are applied to the model. 

STEP 2. A neutral file is created for translation. This creates a preliminary data input file. 
This data file needs to be edited for the non-linear analysis. Certain commands and controls are 
added which are not available in the geometric modelling program. These changes are 
discussed in detail later in the chapter. 

STEP 3. After editing the preliminary data input file, it is ready to be submitted for the 
analysis. (A sample input file is given in appendix A. Appendix B gives a description of 
specific MARC commands.) To submit the job to a supercompter, a job control language (JCL) 
file needs to be prepared. (The JCL file is given in appendix C.) This file contains the 
workspace required and the CPU time required as well as other related details which are 
discussed later. Also a user subroutine file should be present in the same directory as the input 
file. This file suppresses the printing of input in the output file generated after a run. 

After the job has run, four files are created which are as follows: 

1. List file (job.lst) - the output file which gives listing of the results. 

2. Post file (job.post) - the post processing file is used as input to the geometric 
modeling translator. (MARPAT) 

3. Restart File (job.restart) - written if this option is included in the input deck. 
Useful for continuing a run from the last complete run or any complete run for 
which the post data is asked to be written. 
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4. 


Log File (user.O#) - gives details of the job run on the supercompter, giving 
CPU time used and other data. It also gives any SYSTEM errors encountered 
while the job is running. 

If there are any errors the job ends and the list file exists with an EXIT number depending 
on the error. 

STEP 4. The geometric modelling program translator is used to convert the FEA results file 
into readable format. The post file is translated into three files for each increment of the 
analysis which is stored in the post file. These are as follows: 

1. Element File (i#s0.els) - gives the elemental stress values. 

2. Nodal File (i#s0.nod) - gives the nodal stress values. 

3. Displacement File (i#s0.dis) - gives the displacement values. 

4. Log files (marpat.msg and marpat.crd) 

STEP 5. To view the results in the postprocessor, the original neutral file is read to create the 
model. Then results are checked and the files are read as required for elemental, nodal or 
displacement results. 

The flow diagram to show the steps undergone to finally prepare the model for the analysis 
is shown in figure 4.3. A typical MARC input deck is shown in figure 4.4. 

4.5 Changes Made in Preliminary MARC File for Nonlinear Analysis 

Since all of the features available in the MARC program are not available in the geometric 
modeling program, the MARC input file created by the geometric modeling translator must be 
edited to add the following commands and controls for the gear analysis: 

1. SIZING - This value is set to meet the greatest workspace requirement for a given 
problem. The estimate of the greatest workspace requirement is usually based on 
experience. If the SIZING value is not enough for the analysis, the FEA program 
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automatically switches to an out of core solution mode. The out of core solution mode 
is usually not preferred and should be avoided as out of core disk space is required and 
CPU time increases. Initially, when the preliminary FEA input file is obtained SIZING 
has a value of 400,000 words. In this analysis SIZING is changed to 26 Mega words 
(Mw) for the two tooth model and to 28.5 Mw for the seven tooth model. With this 
values the analysis solutions were obtained in-core. Care should be taken in specifying 
this value. This value should be 1.7 Mw less than that specified in the JCL for the 
memory requirement. This accounts for the loading and running of the FEA program. 
This value specified could be 1.7 Mw less than the JCL value, but should never be more 
than the JCL value. 

2. SETNAME - This gives the number of items in defined sets for elements and nodes. 
It usually allots 50 items per set. Since the geometric modeling translator converts each 
named component into a set, a lot of sets are defined with elements and nodes in separate 
sets. This results in excessive name sets being defined. These defined sets are not 
required unless it is necessary to define the nodes and elements, which are difficult to 
identify from the model but are easily defined in these named sets. In that case, only the 
desired name sets are kept and the others are deleted. This not only reduces the input 
file size but also reduces the memory requirements to run the job. Because with an 
increase in name sets, the workspace memory requirement is increased. It is always 
better to reduce the workspace required to permit the job to run within core. Since the 
models used have numerous named components, the preliminary FEA input file has many 
name sets with element and node sets defined separately. The SETNAME has a value 
of 901 initially which is changed to 200 with 50 items per name set i.e., only 4 name 
sets are then kept. This change resulted in the reduction of memory requirements for 
the seven tooth model run. 
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3. DEFINE SETS - These are kept to the minimum and only used if required as is 
discussed in SETNAME change. The only four name sets GEARE, PINIONE, PROJ2E 
and RIME, which define the elements for the respective regions, are used for defining 
individual material properties. 

4. PRINT - This parameter is specified with option 5 to provide output of nodal contact 
information as nodes contact and separate from surfaces during the analysis. 

5. POST - Post processed data is controlled by this option. The type of post process data 
and the rate at which it is recorded is specified. The post codes given are changed to 
17, 131 and 133 which represent VON MISES, MAXIMUM PRINCIPAL and 
MINIMUM PRINCIPAL stresses. The number of increments at which the data is 
recorded is added in column 45. Also 0 in column 20 is changed to 1 to write formatted 
post data. 

6. POINT LOAD - This option, which is given in the beginning of model definition cards, 
is deleted as loads specified in zero (null) increment are ignored. 

7. CONNECTIVITY - In column 15 to suppress element connectivity being printed in the 
output listing a numerical number 1 is added. 

8. COORDINATES - In column 20 to suppress the nodal coordinates being printed in the 
output listing a numerical number 1 is added. 

Additional MARC commands are described in appendix B. 
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4.6 Running Large Finite Element Analysis 

Initially, many trial runs were submitted to study how each of the options behaved. The 
variables in the cards for each option were changed a number of times until an appropriate one 
was found. Considerable effort was spent debugging and optimizing the finite element code 
running on the supercomputer. Because this problem utilized significant computer resources, 
much specialized computer knowledge was needed. This section attempts to document some 
of this experience. 

In the OPTIMIZE option, the method used was changed from 4, the Wave front method 
followed by Grooms, to use option 2, the Cuthill Mckee method. This new method was less 
expensive and saved considerable memory requirements. Also, the ELSTO option has been 
added for this analysis. This reduces the in-core memory requirements below the 28 Mega 
word (Mw) limit defined on the SIZING card. This option was added after a lot of memory 
shortage requirements were faced. Before this problem was resolved more memory was 
required than what could be given in the SIZING card running on the supercomputer. The 
problem finally ran with much less than 200,000 blocks. 

Much experimentation was done on the CONTROL option variables. This option sets limits 
on the number of increments and recycles which may be performed during a nonlinear analysis. 

As discussed earlier, spring elements are used to prevent rigid body modes. The SPRINGS 
option is a list of nodes and their degrees of freedom which are connected by spring elements. 
Nodes were identified on the two bodies where SPRINGS were to be connected. A number of 
trial runs were performed to check the effect of the stiffness value given to the springs. 

The RESTART option is used to save and retrieve analysis data so that the analysis can be 
restarted for additional increments. The variables were adjusted to determine how the loading 
could be continued for the next restart. While performing the initial runs, the restart option was 
used frequently. This was because the number of increments in which the load was applied to 
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reach the final torque in the gears was not appropriate for the solution to converge. However, 
the option was not used in the final analysis, since an efficient the number of load increments 
was eventually found for efficient convergence. 

The loading was specified incrementally with the AUTO LOAD option. AUTO LOAD will 
cause the current load vector to be repeatedly applied (additively) for the number of increments 
requested. The number of increments to be specified in the AUTO LOAD and also the 
CONTROL card sets posed another problem. The actual number of increments to be run were 
specified in the AUTO LOAD card set. The number of increments in the CONTROL card 
were set to the upper limit for the total number of increments allowed for an analysis. This can 
be left blank and let the program take the default maximum allowed. The POINT LOAD 
option is used to apply the load per increment. The FEA program uses incremental values of 
the load which are summed to give the total value that is used in that increment. The last value 
of the incremental load that was input will be used until the new incremental value is read in 
to replace it. 

For contact analysis with AUTO LOAD, a TIME STEP history definition card set must be 
included. Although the time step for each increment is arbitrary for this research analysis, it 
must be included. When rigid bodies are included in a contact analysis, displacements, 
velocities or accelerations are specified and therefore a time step is essential. This problem has 
deformable bodies, not rigid bodies, but the FEA program code still requires a TIME STEP card 
set. An arbitrary time step of 10 units per increment was specified. 

In the CONTACT option, the number of entities present in body 1 and body 2 have to be 
given appropriately depending on the mesh density of the contacting surface. Care has to be 
taken to define the two bodies as body 1 and body 2. Identifying elements present on the 
contacting surface of the bodies becomes complicated for complex models like the spiral bevel 
g ea rs. Many times during this research, the model was optimized in the solid modeler due to 
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some small modifications on the model which caused element numbering to change. This 
necessitated the element numbers to be identified all over again. 

4.7 Convergence Problems 

For a nonlinear analysis, many load steps or increments must be used to reach the final or 
fully loaded state. In the process, one can encounter some convergence problems. 

A torque of 9508 in. lbs. was applied to the gear. This torque was input as a series of point 
loads. Since this is a nonlinear problem, the entire torque should not be applied in one 
increment. Several analysis (about 60) were run to determine the sensitivity of the structure 
meeting the convergence criteria as a function of the size of the applied load. 

The FEA program solves the non-linear problem on an incremental basis and iterates within 
each increment until the convergence criteria are met. These criteria are specified in the 
CONTROL data block. Too many iterations within an increment, is an indication that too 
much loading is being attempted for that load step. Initially the gears are separated by a small 
gap. At this point the total structure has very little stiffness. Applying any torque causes the 
gear to rotate and contact the adjoining pinion. Very small load increments should be applied 
otherwise the matrix updated for this soft structure will not converge. After the gear contacts 
the stiffness of the structure increases drastically. The amount of loading can then be increased 
until the required amount of torque is achieved. For the seven tooth model, a total of 8 
increments were used to reach the desired torque, while for the two tooth model a total of 10 
increments were used. 

While performing the trial analysis some increments kept repeating as the convergence 
criteria was not met. Also it was observed that load was automatically decreased by the 
program. The conclusion was drawn that the load applied was too high for the solution to 
converge, which was seemed correct as the program itself was decreasing the load values. This 
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gave an indication that the load step needed to be reduced. 

To determine the correct load to be applied many trial analysis were performed. The load 
was reduced from about 100 lbs. to 5 lbs. but the solution was still not converging. Each time 
different nodes kept coming into contact and disturbed the convergence requirement. The nodes 
in contact were found to lie on the border of the defined contact region. Since nodes other than 
those defined as contact region nodes were trying to contact, convergence was affected. 

The next step taken was to define more surface nodes on the contact region (CONTACT 
CARD). After this the convergence criteria was changed from displacement to force residuals 
(CONTROL CARD). This eventually led to the convergence of the analysis. 

A summary table of the computer runs and the incremental loads given to the seven tooth 
model and two tooth model is given in Table III. 


TABLE III: SUMMARY OF LOADING AND CPU TIME 


SEVEN TOOTH MODEL 

Increment Number 

Total CPU time (sec) 

Load increment (lbs.) 

Total Load (lbs.) 

1 

437.5 

12 

12 

2 

1403.37 

2 

14 

3 

1968.24 

785 

799 

4 

2541.77 

785 

1584 

5 

3115.01 

785 

2369 

6 

4287.06 

785 

3939 

7 

4886.40 

785 

4274 
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TWO TOOTH MODEL 

Increment Number 

Total CPU time (sec) 

Load Increment (lbs.) 

Total Load (lbs.) 

1 

121.11 

2 

2 

2 

150.38 

2 

4 

3 

181.13 

2 

6 

4 

212.98 

2 

8 

5 

309.7 

2 

10 

6 

405.28 

2 

12 

7 

775.31 

845 

857 

8 

856.49 

845 

1702 

9 

939.53 

845 

2547 

10 

1097.90 

845 

3392 


4.8 Supercomputer Requirements 

Supercomputers have various queues with different memory and CPU time requirements as 
shown in Table IV. Each of the queues has either 1 or 2 jobs running. A small job with 
SIZING of around 400,000 words usually takes less than 300 secs of CPU time and can be in 
the smallest queue called debug. But as the size of a problem increases the CPU requirements 
increase which necessitates a job to be submitted in a higher queue. The 2 jobs analyzed in this 
research required 26 Mw and 28.5 Mw of memory space and about 20 - 80 minutes of CPU 
time. The turnaround from each job was very slow. A job on an average would take a day 
and a half in the queue before running. After a lot of experimentation, the job was first 
submitted in the 300 second slot. Because of this, errors if any would show up in 10 minutes 
rather than in 2 days. 
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TABLE IV: TYPICAL INPUT QUEUES FOR SUPERCOMPUTER 




1 Nisht Queue Limits 

Weekend Queue Limits 


Queue 

Name 

warn 


Max Jobs Running 

Max Jobs Running 

Max Jobs Running 

■ 

ESsai 

Class User 

Class User 

Class User 

debus 

300 

30.2 

3 

i 

2 

1 

1 

1 

2 

am 

1.200 

■i 

2 

i 

2 

1 

1 

1 

2 

warn 

3.600 

4.2 

2 

i 

2 

1 

1 

1 

2 

!■ 

7,200 

4.2 

2 

i 

2 

1 

1 

1 

1 


3.600 

8.2 

2 

i 

2 

1 

1 

1 

2 

a 7 

7.200 

8.2 

2 

i 

2 

1 

1 

1 

1 

a 8 

10.800 

8.2 

2 

i 

1 

1 

1 

1 

1 

q 9 

7.200 

16.2 

i 

i 

1 

1 

1 

1 

1 

a’ 0 

14,400 

16.2 

i 

i 

1 

1 

1 

1 

1 

a" 

21.600 

16.2 

i 

i 

1 

1 

I 

1 

1 

a 12 

10.800 

30.2 

i 

i 

1 

1 

1 

1 

1 

a 1J 

21.600 

30.2 

. 


1 

1 

1 

1 

1 

— 

21.600 

60.2 

- 

- 

- 

- 

1 

1 

1 


Day Queues: 0800-1700 Mon. thru Thurs. & 0800-2200 Fri. 
Night Queues: 1700-0800 Mon. thru Fri. 

Weekend Queues: 2200 Fri. thru 0800 Mon. 

Allow 0.7 Mw for System Memory Overhead 


As mentioned earlier, the space requested in the job control data (JCL) set should be 1.7 Mw 
higher than that given in the SIZING card. This is required for the FEA program to be loaded. 
A detailed listing of the JCL is given in appendix C. 

If the job is required to be restarted, the input and output RESTART files should be 
mentioned in the correct format. This option was utilized in this research since the turnaround 
was very slow. Only a few increments were requested for the analysis and were saved in the 
restart file. By reading the restart file, the job could be restarted from the previous increment. 

Computing the amount of workspace required by the FEA program is a complex function 
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of many variables. The most efficient method is to select a large enough workspace to handle 
a variety of runs, without sacrificing efficiency of wasting core space. 

Both in-core and out-of-core data storage schemes are available in the FEA program. 
Elements may also be stored out-of-core if the ELSTO option is used. The FEA program 
chooses the solution type automatically. 

The in-core solution technique is used when the workspace size specified in the SIZING card 
is larger than the total workspace needed in the in-core matrix. When the workspace required 
is too large, program uses the out-of-core solution techniques and show how much space each 
nodal row requires, the number of nodal rows per buffer and how large each auxiliary file would 
be. The buffer size can be increased only by changing allocation on the SIZING card. The 
amount of size to be given is based on experience. For large problems, the exact workspace 
requirements can be determined without actually running the job by inserting a STOP parameter 
card after the workspace is allocated. 

The frequency for writing restarts and POST data should be low to avoid disk space 
problems since these files are very large. Initially these files were written after every 10 
iterations. Finally they were written for every 2 iterations when the number of increments to 
apply the load was reduced for the job to run. 
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CHAPTER V 


RESULTS 

Various jobs were successfully run to predict the contact path and the load sharing for the 
double tooth contact region as the gears roll through mesh. 

5.1 Elliptical Stress Distribution 

The stress distribution in the contact region was found to be elliptical. Figure 5. 1 shows a 
typical elliptical contour obtained using elemental stress values on the pinion surface. Figure 

5.2 shows typical pinion contact stresses with nodal values. Because of the large nodal stress 
gradient, the nodal stress ellipse is slightly distorted when compared to the ellipse contour with 
elemental stresses. Note that only the pinion contours and the stress values are discussed in this 
research, since the gear teeth share a similar load distribution. Also note that only the 
minimum principal stresses are recorded at the contact region for the study. 

5.2 Gap Element and Automated Contact Analysis Comparison 

Contact stresses on spiral bevel gears were studied by researchers with gap elements in 
references 16 and 17 on a similar seven tooth model. Comparison of the results with automated 
contact analysis will now be presented. Both models contained the same mesh density, 
boundary conditions, material properties and loading. The nodal stress results of pinion tooth 
#1 obtained from gap elements and automated contact analysis are as shown in figures 5.3 and 
5.4. Note that both the contours show the highest concentrated stress value at the same node. 
A comparison of the results of these two runs are as follows: 
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GaD Element 

Automated Contact 
Analysis 

Max Nodal stress 

-296,410 (psi) 

-291503 (psi) 

CPU time (approx) 

30 min 

80 min 

Elemental stress 

-84,761 (psi) 

-113,577 (psi) 

Gap elements closed or 
(nodal points with contact) 

4 

8 

No. of increments 

4 

8 


The number of contact nodes at the contact region in the automated contact FEA program was 
higher than identified by the gap element FEA program. With the pinion considered body 1 
(contactor body) and the gear considered body 2 (contacted body), eight nodes contacted as 
shown above. With body 1 and body 2 switched, 16 nodes were found to have contact. 
Presumably this sort of discrepancy occurred because the mesh was too coarse. 

5.3 Seven Tooth Model Results 

As discussed earlier, the seven tooth model pinion was rotated from +6 degrees to -30 
degrees about the X axis and the gear was rotated from +2 to -12 degrees about the Z axis. 
As these gears were rolled there was a shift in the contact region. It was observed that as one 
pinion tooth goes out of mesh and the load on it reduces, the other pinion tooth shares more load 
and as it starts coming into contact. The elliptical stress contours for pinion tooth #1 and pinion 
tooth #2 while they are rotated in mesh with the respective gear teeth (for all the rotations stated 
above) are shown in figure 5.5. 

Figures 5.6 to 5. 19 show plots of nodal stresses for pinion tooth #1 and #2 for all positions. 
Shift in the contact nodes while the gearset rolls through mesh are shown in figures 5.20 and 
5.21. These nodes are obtained from the output file. The contact forces on the pinion and the 
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gear nodes which are listed at the end of the last increment are studied to identify the nodes in 
contact. The contact node density, or all the nodes that contact during meshing, is as shown 
in figure 5.22. 

5.4 Two Tooth Model Results 

The two tooth model was generated with a finer mesh with 24 nodes along the length of the 
tooth and 12 along the height of the tooth. 

As discussed earlier, the two tooth pinion is rotated in six degree increments to the following 
positions: +6, 0, and -6 degrees to observe the shift in the contact ellipse. The direction of 
rotations were similar to that of the seven tooth model. The element and nodal stress results 
for the model rotated by -6 degrees are shown in figures 5.23 and 5.24. In this model due to 
the flexible hub, the gear started sliding over the pinion by a large amount resulting in more 
nodes in contact. For this reason the contact ellipse was sliding towards the edge of the pinion 
and the ellipse contour became distorted. A summary of the seven tooth and two tooth models 
with reference to modeling and analysis are given in Table V. 


TABLE V: SUMMARY OF FEA ANALYSIS 


FEATURES 

7 - TOOTH MODEL 

2 - TOOTH MODEL 

No. of elements 

8793 

3116 

No. of nodes 

11261 

4452 

No. of degrees of freedom 

33748 

13321 

"In-core" Memory required (Mw) 

28.5 

26 

CPU time (seconds) 

4886.40 

1092.90 

No. of increments 

8 

10 

Contact tolerance (in.) 

0.0002 

0.0002 

Torque applied lbs. in. 

9508 

9508 
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5.5 Contact Path 

The contact path is defined as the path of the point of maximum force on the gear tooth 
while it rolls in and out of mesh. These values are determined for the pinion using the output 
listing which gives the contact force in the final increment before converging. Referring to 
Figure 5.5, it should be noted that while the contact ellipse in pinion tooth #2 shifts from the 
toe of the tooth to the center of the tooth, the contact ellipse in pinion tooth #1 shifts from the 
center of the tooth to the heel of the tooth. The contact path for pinion tooth #1 and pinion tooth 
#2 are plotted in figure 5.25. The two curves of pinion tooth #1 and pinion tooth #2 do not 
overlap because the mesh density is too coarse and also the rotations are not small enough. 

The force at a particular node is taken to be the square root of the sum of the squares of 
the forces in the X,Y and Z directions. The maximum force is obtained for all the rotations 
of the model and are plotted on the pinion. The two tooth model showed a shift in the contact 
path due to flexibility effects. The hub region for this model was very flexible and had 
deformed a lot before the gear pair came into contact. 

5.6 Load Sharing in Spiral Bevel Gears 

The load sharing was analyzed in the double tooth contact region of the seven tooth 
model. The total forces on pinion tooth #1 and pinion tooth # 2 at each of the rotations from +6 
degrees to -30 degrees for the seven tooth model are calculated. This is done using the contact 
table obtained from the output listing of any run. The total contact force on each pinion is 
determined by 

/(ZF,)*’ ♦ (ZF y ) 2 . < 14 > 

Table VI shows the load on pinion tooth #1 and pinion tooth # 2 calculated from various runs 
while the gears were rotated. The load on the tooth as a function of rotation is also plotted as 
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shown in figure 5.26 


TABLE VI: MAGNITUDE OF CONTACT FORCES 
IN PINION 1 AND PINION2 ACROSS THE CONTACT REGION 


Degree of Rotation 

PINION TOOTH #1 

PINION TOOTH n 

of Pinion 

F x 

HU 

F z 

F t , 

F x 

H 

F z 

F* 

+6 

2870 

3294.4 

473.6 

4397.4 

0 

0 

0 

0 

0 

2710.4 

2759.0 

447.0 

3893.0 

154.8 

164.9 

29.9 

226 

-6 

2323.7 

2032.0 

391.17 

3110.98 

513.0 

515.0 

83.5 

767.7 

-12 

1655.0 

1908.4 

295.3 

2542.9 

1008.9 

1093.4 

164.5 

1495.91 

-18 

1114.0 

1267.1 

240.5 

1704.07 

1582.2 

1655.1 

175.5 

2296.16 

-24 

646.3 

763.2 

135.7 

848.2 

1873.3 

2019.4 

343.68 

2775.2 

-30 

200.91 

252.6 

45.8 

324.8 

2156.0 

2326.0 

450.4 

3203.33 
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5.7 Fillet Stresses 


The fillet stresses are plotted for the seven tooth model with load sharing. A typical 
contour of these stresses for one particular rotation is shown in figure 5.28. Table VII gives the 
maximum principal nodal stresses at various nodes and at various roll angles. The location of 
these nodes are given in figure 5.27. These values are plotted as a function of roll angle as 
shown in figure 5.29. 


TABLE VII: FILLET STRESSES AS A FUNCTION OF ROTATION 
FOR DIFFERENT NODES 


PINION 1 

|j 

Nodes- > 

883 

871 

856 

1745 

1733 

1718 1 

Rotations | 

+6 

75966.13 

138493.3 

32738.5 

39224 

18491 

15248 

0 

44804.1 

109165.1 

37978.5 

45357.8 

16021.8 

12918.5 

-6 

24030.4 

66600.7 

29362.4 

67333.4 

17878.6 

10812.1 

-12 

10272.5 

21690.6 

12490.4 

74458.8 

28887.1 

9505.6 

-18 

3721.5 

1312.2 

5926.3 

81997 

55723 

11607.1 

-24 

414.4 

15953.0 

13148.4 

90294.8 

67987.5 

14639.2 

-30 

1806.3 

29486.5 

16799 

37915.3 

81093.5 

27288.9 
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CHAPTER VI 


SUMMARY AND CONCLUSIONS 

A procedure for predicting the contact path, load sharing, contact stresses and fillet stresses 
of spiral bevel gears is presented. The method incorporates the following steps: 

1. A model was developed using the kinematics of the manufacturing process, the machine 
settings, and the design data for the gearset. The model was generated in a geometric 
modeling program. 

2. Automated contact analysis option in a nonlinear finite element analysis program was 
used to perform the analysis. 

3. Two different models were analyzed. Model I consisted of three gear teeth and four 
pinion teeth in mesh. The nodal mesh density on the tooth surface was 10 by 10. 
Model II consisted of one pinion and one gear tooth in mesh with a nodal tooth surface 
mesh density of 24 by 12. 

4. These models were analyzed by rotating the gears in angular increments to predict the 
contact path. Model I was rotated a total of 36 degrees and Model II was rotated a total 
of 12 degrees. 

5. The load sharing was determined from the contact forces on the nodes at the end of the 
converged solution for each rotational position. 

The initial FEA results for Model I at the contact region compared favorably with the results 
by earlier researchers using gap elements. It was observed that in the double contact zone, 
when the contact ellipse in the (i)th pinion shifted from the center of the gear tooth to the heel 
of the gear tooth, the contact ellipse in the (i+ l)th pinion shifted from the toe of the gear tooth 
to the center. The load for each tooth was calculated as the sum of the contact forces on 
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various nodes in the contact region. The number of nodes in contact changes as the gear rolls 
through mesh. 

There was a large stress gradient between adjacent nodes in the contact region. This 
indicates a need for mesh refinement. 
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Figure 2.2 Gear and Pinion Orientations to Mesh with Each Other 
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Figure 2 . 3 


3-D Model of the Spiral Bevel Gears 
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Figure 3.1 The Full Newton-Raphson (N-R) Method 
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Figure 3.2 The Modified Newton-Raphson Method 
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BODY 1 


Contact Tolerance 



When contact is detected, Multi-Point 
Displacement Constraint Ties are 
introduced to prevent body penetration 


BODY 2 



Figure 3.4 Contact Algorithm 
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Figure 4.2 Two tooth model with boundary conditions 
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Figure 4.3 Flow diagram for the MARC analysis 
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Figure 4.4 MARC Input Deck 
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PSI 



A -6664. 
B -18543. 
C -30422. 
D -42302. 
E -54181. 
F -66060. 
G -77939. 
H -89819. 
I -101698. 
J -113577. 
K -125457. 
L -137336. 
M -149215. 
N -161094. 
O -172974. 


Figure 5.1 A typical elliptical contour in the pinion 

tooth with elemental stresses in seven tooth 
model (PSI) 
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PSI 



A 85251. 
B 22459. 
C -40334. 
D -103126. 
E -165918. 
F -228711. 
G -291503. 
H -354295. 


Figure 5.2 A typical elliptical contour in the pinion tooth 
with nodal stresses in seven tooth model 
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PSI 



-17629. 

-57455. 

-97281. 

-137107. 

-176933. 

-216759. 

-256584. 

-296410. 

-336236. 

-376062. 

-415888. 

-455714. 

-495540. 

-535366. 

-575191. 


Figure 5.3 Nodal Stress Result on Pinion Obtained From the Gap 
Element Solution and Seven Tooth Model 
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PSI 


A 

B 

C 

D 

E 

P 

G 


85251. 

22459. 

-40334. 

-103126. 

-165918. 

-228711. 

-291503. 



Figure 5.4 Nodal stress result on pinion obtained from 

automated contact analysis in seven tooth model 
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Figure 5.5 Elemental Stresses for Pinion Tooth #1 and Pinion 
Tooth #2 as They Roll In and Out of Mesh for 
Rotations from +6 to -30 Degrees in Seven Tooth Model. 
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Pinion Tooth #1 


-12 Degree rotation Pinion Tooth #2 





Figure 5.5 (continued) 
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PSI 


A 109199. 



Figure 5.6 Nodal stresses on pinion tooth #1 after it is 
rotated by +6 degrees in the seven tooth model 
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PSI 



A 22611. 
B 16234. 
C 9857. 
D 3480. 
E -2897. 
F -9274 . 
G -15651. 
H -22028. 
I -28405. 
J -34782. 
K -41159. 
L -47536. 
M -53913. 
N -60290. 
O . -66666 . 


Figure 5.7 Nodal stresses on pinion tooth #2 after it is 
rotated by +6 degrees in the seven tooth model 
(No Contact) 
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PSI 



A 85251. 
B 22459. 
C -40334. 
D -103126. 
E -165918. 
F -228711. 
G -291503. 


Figure 5.8 Nodal stresses on pinion tooth f 1 after it is 
rotated by 0 degrees in the seven tooth model 
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PSI 



A 94082. 
B 50628. 

C 71 74. 

D -36280. 
E -79734. 
F -123188. 
G -166642. 
H -210096. 
I -253550. 
J -297004. 
K -340458. 
L -383912. 
M -427366. 
N -470820. 
O -514273. 


Figure 5.9 Nodal Stress on Pinion Tooth #2 after it is rotated 
0 Degrees in the Seven Tooth Model 
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PSI 



A 101003. 
B 27684. 
C -45634. 
D -118952. 


Figure 5.10 Nodal stresses on pinion tooth #1 after it is 
rotated by -6 degrees in the seven tooth model 
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PSI 



46780 . 

11692 . 

- 23396 . 

- 58484 . 

- 93572 . 

- 128660 . 

- 163748 . 

- 198836 . 

- 233924 . 

- 269012 . 

- 304100 . 

- 339188 . 

- 374276 . 

- 409364 . 

- 444452 . 


Figure 5.11 Nodal Stress on Pinion Tooth #2 after it is rotated 
-6 Degrees in the Seven Tooth Model 
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PSI 



A 84411. 
B 17923. 
C -48566. 
D -115054. 


Figure 5.12 Nodal stresses on pinion tooth #1 after it is 

rotated by -12 degrees in the seven tooth model 
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PSI 



A 29291. 

B 11202. 

C *6887. 

D -24976. 

E -43066. 

F -61155. 

G -79244. 

H -97334. 




j -133512. 
K -151601. 
-169691. 

Ij 

„ -187780. 

M 

„ -205869. 

N 

Q -223959. 


Figure 5.13 Nodal Stress on Pinion Tooth #2 after it is rotated 
-12 Degrees in the Seven Tooth Model 
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PSI 



A 63612. 
B 2163. 
C -59286. 
D -120735. 


Figure 5.14 Nodal stresses on pinion tooth #1 after it is 

rotated by -18 degrees in the seven tooth model 
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PSI 



59076 . 

31539 . 

4003 . 

- 23534 . 

- 51071 . 

- 78607 . 

- 106144 . 

- 133681 . 

- 161217 . 

- 188754 . 

- 216291 . 

- 243828 . 

- 271364 . 

- 298901 . 

• 326438 . 


Figure 5.15 Nodal Stress on Pinion Tooth #2 after it is rotated 
-18 Degrees in the Seven Tooth Model 
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PSI 


A 41824. 
B -20362. 




C -82549. 
D -144735. 


Figure 5.16 


Nodal stresses on pinion tooth #1 after it is 
rotated by -24 degrees in the seven tooth model 
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PSI 



A 19971. 

B -2764. 

c -25499. 

-48234. 

E -70968. 

F -93703. 

G -116438. 

H -139173. 

x -161908. 

j -184643. 

K -207378. 

-230112. 

L 

M -252847. 
„ -275582. 

O -298317. 


Figure 5.17 Nodal Stress on Pinion Tooth #2 after it is rotated 
-24 Degrees in the Seven Tooth Model 
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PSI 


A 21750. 
B -40817. 
C-103384. 



Figure 5.18 Nodal stresses on pinion tooth #1 after it is 

rotated by -30 degrees in the seven toot h m odel 
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PSI 


A 24987. 

B 10637. 

c -3713. 

D -18063. 

E -32413. 

F -46763. 

G -61113. 

H -75463. 

X -89813. 

J -104163. 

K -118513. 

L -132863. 

M -147213. 

161563. 


175913. 








mtacting Nodes 


St 
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Coi 


Nodal Contact Points on Gear. Shift xn the 
contact nodes on gear teeth at selected 
rotations during meshing for the seven tooth 
model 







Figure 5.22 Contact Node Density for the Seven Tooth Model 
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PSI 



A -4834. 
B -16283. 
C -27731. 
D -39180. 
E -50629. 
F -62077. 
G -73526. 
H -84975. 
I -96424. 
J -107872. 
K -119321. 
L -130770. 
M -142219. 
N -153667. 
O -165116. 


Figure 5.23 Elemental stress results for two tooth model 
at -6 degree rotation 
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PSI 



A 198082. 
B 128166. 
C 58251. 
D -11665. 
E -81580. 
F -151496. 
G -221411. 


Figure 5.24 Nodal stress results for two tooth model at 
-6 degree rotation 
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Figure 5.25 Contact Path of Seven Tooth Model in Pinion Tooth #1 
and Pinion Tooth #2. Shown is the Location of the 
Highest Contact Force For Each Rotation Indicated. 
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FORCE ON PINION TOOTH IN LBS, 


LOAD AS A FUNCTION OF ROTATION 



.80 -60 -40 -20 0 20 40 60 

ROLL ANGLE IN DEGREES 

Figure 5.26 Load on Tooth as a Function of Rotation 


Figure 5.26 Load on Tooth as a Function of Rotation 
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Figure 5.27 Location of Nodes Identified in Table VII. 
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Pinion Tooth #1 




Figure 5.28 Typical fillet stress contours. 

Shown for -6 degree rotation of pinion 
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APPENDIX A 


MARC INPUT DATA 


TITLE BULK DATA CARDS PRODUCED BY PATMAR RELEASE 3.1A 
TITLE 20-Jul-93 15:14 

TITLE OPTIMIZED 7 TEETH MODEL 
S *** PARAMETER CARDS *** 

$.. . ! 1 ! 2 ! 3 ! 4 ! 5 ! .. . 


.8 


SIZING 28500000 ... this allocates workspace storage memory 

SETNAME 200 ... defines max. no of items in DEFINE sets 

ELEMENTS 7 ... element type is 8 noded hexahedral 

PRINT, 5 — gives extra information about contact nodes 

ELSTO . . . reduces in-core memory by storing elements out-of-core 

END . . . signifies end of parameter cards 

S *** CONTROL CARDS FOR POST-PROCESSOR TAPE *** 

S SPECIFY POST CODES AND LABELS FOR EACH ELEMENT VARIABLE. 

$...! 1 ! 2 ! 3 ! 4 ! 5 ! 6 ! 7 ! 8 

POST 

$ 1 in 20th column gives formatted pose file which is readable on different m/c's 
S 4 shows the frequency at which POST data is written 


3 11 4 

17 VON MISES STRESS 

131 MAXIMUM PRINCIPAL STRESS 

133 MINIMUM PRINCIPAL STRESS 

S *** MODEL DEFINITION CARDS *** 

$ SPECIFY DISPLACEMENTS AND THE ASSOCIATED DOFS AND NODES FOR EACH SET. 



1 . 

t 

2. 

1 

3. 

t 

, . . . 4 

! 5 ! 6 ! 7 . . ., 

FIXED 

DISP 








3 

0.0 







three sets of displacement data 
prescribed displacement 

3 







. . • 

Z direction constrained 

2952 

2961 

2992 

3016 

8500 

8581 

9680 

9761 

. . . list of nodes on gear 


0.0 


0.0 


0.0 



1 

2 

3 





. . . 

X/Y/Z directions constrained 

7872 

0.0 


0.0 


0.0 



pivot node on hub 

1 

2 

3 







259 

263 

364 

380 

2317 

2333 

2528 

2561 

. . . list of pinion nodes 

$ — 

A NAMED SET MUST BE 

DEFINED BEFORE ITS 

NAME IS USED IN THE INPUT DECK. 

$...!. 

1. 

! 

2. 

t 

3. 

t 

4 ! 



SElement set GEARE is th list of elements of named 

component 

GEAR 

in PATRAN 


DEFINE ELEMENT SET 


GEARE 









1837 1838 1839 1840 1841 

1842 

1843 1844 

1845 

1846 

1847 

1848 

1849 

1850 

1851 

C 

7426 7427 7428 7429 7430 
7441 7442 7443 

7431 

7432 7433 

7434 

7435 

7436 

7437 

7438 

7439 

7440 

c 

$De£ines PINION elements 
DEFINE ELEMENT SET 


PINIONE 









1 2 3 4 5 

6 

7 8 

9 

10 

11 

12 

13 

14 

15 

c 

1831 1832 1833 1834 1835 
DEFINE ELEMENT SET 

1836 

PROJ2E 









5023 5024 5025 5026 5027 

5028 

5029 5030 

5031 

5032 

5033 

5050 

5051 

5052 

5053 

c 

6679 6680 6681 6682 6683 

6684 

6685 6686 

6687 








DEFINE ELEMENT SET 


RIME 









2116 2117 2118 2119 2120 

2121 

2122 2123 

2124 

2125 

2126 

2127 

2128 

2129 

2130 

c 

8779 8780 8781 8782 8783 

8784 

8785 8786 

8787 

8788 

8789 

8790 

8791 

8792 

8793 


S SPECIFY CONNECTIVITY FOR 

EACH ELEMENT. 








S__ .! 1 ! 2 ! . 

3. 

i 4. 

t 

5. 

» 

6. 

1 

7. 

I 

. _ a 

CONNECTIVITY 


Sshows total no. of elements in the model, 1 suppresses connectivity printout 
8793 1 
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3 

7 

54 

53 

3 

4 

59 

58 

8 

9 

4 

7 

55 

54 

4 

5 

60 

59 

9 

10 

5 

7 

57 

56 

6 

7 

62 

61 

11 

12 


8791 71124611248112581125211247112491126011256 

8792 711251 1125711254 112501 12 5S1126111259112S3 

8793 71125211258112571125111256112601126111255 

$ SPECIFY COORDINATES FOR EACH NODE. 

$...! 1 ! 2 ! 3 ! 4 ! 5 ! 6 ! 7 ! 8 

COORDINATES 

$# of directions per node, total no. of nodes, suppression of coordinate listing 
611261 1 

1 -3.429364 -0.538929 1.299387. 

2 -3.429364 -0.553653 1.292881 

3 -3.429364 -0.568376 1.286376 

4 -3.429364 -0.S83099 1.279870 

5 -3.429364 -0.597823 1.273365 

11261 -1.094289 -0.899709 1.786157 

Sdefines material properties 
ISOTROPIC 

3 ... three sets of material , properties 

1, VON MISES 

3.0000E+07,0.3 ... defines Young's modulus, poisson ratio 

Slist of elements not nodes . . . 

SLIST OF GEAR AND PINION ELEMENTS . . . 

GEARE AND PINIONE — elements defined by ELEMENT SET names 

2, VON MISES 

3. 000OE+10, 0 . 3 
5RIM ELEMENTS... 

RIME 

3, VON MISES 
6.0000E+13, 0.3 

SPRJECTION IN FRONT OF THE HUB.. 

PROJ2E 

SPRINGS 

258.1.2671.1.100 ... 1st node, dof, 2nd node, dof, stiffnes 

258.2.2671.2.100 

258.3.2671.3.100 

364.1.2626.1.100 

364.2.2626.2.100 

364.3.2626.3.100 

2528.1.5512.1.100 

2528.2.5512.2.100 

2528.3. 5512. 3. 100 

2317.1.5541.1.100 

2317.2.5541.2.100 

2317.3.5541.3.100 
CONTACT 

2,450,4SO no. of contacting bodies and entities on their surface 

,0.0002,, contact tolerance defined 

1,0 definition of 1st body 

t t t 

o:, 

Slist of gear elements not nodes on the contacting surface....... 

$GE_4 BACK FACE ELEMENTS.. 

2902 TO 2904 BY 1 AND 3109 TO 3111 BY 1 AND 3316 TO 3318 BY 4 AND C 

3523 TO 3525 BY 1 AND 3730 TO 3732 BY 1 AND C 

2666,2670, AND 2673 TO 2685 BY 4 AND 2686 TO 2688 BY 1 AND C 

2881 TO 2901 BY 4 AND 3088 TO 3108 BY 4 AND C 

3295 TO 3315 BY 4 AND 3S02 TO 3522 BY 4 AND C 

3709 TO 3729 BY 4 AND 3745 TO 3771 BY 1 AND C 
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$GE_S BACK FACE ELEMENTS . . . 

4X80 TO 4182 BY 1 AND 4450 TO 4452 BY 1 AND C 

4702 TO 4704 BY 1 AND 4945 TO 4947 BY 1 AND C 

5215 TO S217 BY 1 AND 5485 TO 5487 BY 1 AND C 

4159 TO 4179 BY 4 AND 4429 TO 4449 BY 4 AND C 

4681 TO 4701 BY 4 AND 4924 TO 4944 BY 4 AND C 

5194 TO 5214 BY 4 AND 5464 TO 5484 BY 4 AND C 

5734 TO 5739 BY 1 AND 5806 TO 5811 BY 1 AND C 

5740 TO 5742 BY 1 AND 5815 TO 5820 BY 1 

2/0 ... defi.nci.on of 2nd body 

t t t 
0., 

Slist of pinion elements not nodes on the contacting surface 

$PI_1 FRONT FACE ELEMENTS . . 

4 TO 16 BY 4 AND 40 TO 52 BY 4 AND 76 TO 88 BY 4 AND C 
112 TO 124 BY 4 AND 148 TO 160 BY 4 AND 184, 195, 199, 203, C 
282, 291, 300, 303, 362, 369, 376, 383, 460, 465, 470, 475, C 
20,27,32,35,36,56,63,68,71,72,92,99,104,107,108, AND C 
128, 135, 140, 143, 144, 164, 171, 176, 179, 180, AND C 

207, 214, 219, 222, 223, 306, 309, 314, 317, 318, AND C 
385,387,389,392,393, AND 480 TO 484 BY 1 AND C 
SPI_2 FRONT FACE ELEMENTS . . . 

508 TO 520 BY 4 AND 580 TO 592 BY 4 AND 652 TO 664 BY 4 AND C 
726 TO 738 BY 4 AND 796, 807, 811 , 815, 903, 912, 915, 988, 995, C 
443,450,455,458,459,524,531,536,539,540, AND C 
596,603,608,611,612,668,675,680,683,684, AND C 
742,749,754,757,758,819,826,831,834,835, AND C 
918,921,928,929,930,997,999,1001,1004,1005, AND C 
1092 TO 1096 BY 1 
CONTACT TABLE 

1 ...one set of bodies defined 

1 . . .bodyl is contacting 

2 . . .body 2 is contacted 

OPTXMX2E, 2, , , , 1 . . . Cuthill Mckeee bandwidth optimizer is used 

20 , 

CONTROL ...sets controls for iterations 

SI of increments, max. I of recycles, min. I of recycles, residual checking invoked,.. 
, 5, 0, 1, 0, , 

.15 

RESTART ...optional restart command 

1,2 

PRINT ELEM 
1,2 

STRAIN, STRESS 
1 TO 8793 

PRINT NODE 

1,2 

LOAD, REAC, TOTA, STRESS 
SCONSTRAINT POINTS . . . 

7872,2952,2961,2992,3016,8500,8581,9680,9761, C 
259,263,364,380,2317,2333,2528,2561, AND C 
SPI_1 FRONT FACE NODES . . . 

5 TO 20 BY S AND 55 TO 70 BY 5 AND 105 TO 120 BY 5 AND C 

155 TO 170 BY 5 AND 205 TO 220 BY 5 AND 255, 268, 273, 278, C 
384, 395, 406, 410, 493, 502, 511, 520, 632, 639, 646, 653, 751, AND C 
25 TO 225 BY 50 AND 34 TO 234 BY 50 AND C 

41 TO 241 BY SO AND 46 TO 246 BY SO AND C 

49 TO 249 BY 50 AND 50 TO 250 BY SO AND C 

283, 292, 299, 304, 307, 308, 414, 418, 425, 430, 433, 434, 761, 756, C 
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523,526,529,534,537,538, AND 660 TO 670 BY 2 AND C 
770 TO 780 BY 2 AND C 
$PI_2 FRONT FACE NODES... 

585 TO 605 BY 5 AND 614,621,626,629,630, AND C 

705 TO 725 BY 5 AND C 

805 TO 825 BY 5 AND C 

90S TO 925 BY 5 AND C 

734 TO 934 BY 100 AND C 

741 TO 941 BY 100 AND C 

746 TO 946 BY 100 AND C 

749 TO 949 BY 100 AND C 

750 TO 950 BY 100 AND C 
1009 TO 1029 BY 5 AND C 

1038, 104S, 1050, 1053, 1054, 1117, AND C 

1130 TO 1145 BY 5, AND 1154,1161,1166,1169,1170, AND C 

1246, 12S7, AND 1268 TO 1280 BY 4 AND 1287,1292,1295,1296, AND C 

13SS, 1364, 1373, 1382, 138S, 1388, 1391, 1396, 1399, 1400, AND C 

1494, 1501, 1508, 1515, AND 1522 TO 1528 BY 2 AND 1531,1532, AND C 

1613, 1618, AND 1623 TO 1638 BY 5 AND 1639 TO 1642 BY 1 AND C 

$GE_4 BACK FACE NODES... 

3678 TO 3681 BY 1 AND 3948 TO 3951 BY 1 AND C 

4198 TO 4201 BY 1 AND 4448 TO 4451 BY 1 AND C 

4698 TO 4701 BY 1 AND 4948 TO 4951 BY 1 AND C 

3652 TO 3677 BY 5 AND 3922 TO 3947 BY 5 AND C 

4172 TO 4197 BY 5 AND 4422 TO 4447 BY 5 AND C 

4672 TO 4697 BY 5 AND 4922 TO 4947 BY 5 AND C 

4972 TO 5011 BY 1 AND C 

$GE_S BACK FACE NODES... 

5488 TO 5491 BY 1 AND 5828 TO 5831 BY 1 AND C 

6148 TO 6151 BY 1 AND 6448 TO 6451 BY 1 AND C 

6778 TO 6781 BY 1 AND 7108 TO 7111 BY 1 AND C 

5462 TO 5487 BY 5 AND 5802 TO 5827 BY 5 AND C 

6122 TO 6147 BY 5 AND 6422 TO 6447 BY 5 AND C 

6752 TO 6777 BY 5 AND 7082 TO 7107 BY 5 AND C 

7412 TO 7421 BY 1 AND 7512 TO 7541 BY 1 

ERROR ESTIMATE 

1,1 

END OPTION ....terminates model definition cards 

SHistory or Load Definition Cards 

POINT LOAD 

0 . 0 , 12 . 0 , 0.0 

5141 

AUTO LOAD 
1 

TIME STEP 

10.0 

CONTINUE 
POINT LOAD 

0 . 0 , 2 . 0 , 0.0 

5141 

AUTO LOAD 
1 

TIME STEP 

10.0 

CONTINUE 
POINT LOAD 

0.0,785.0,0.0 
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5141 

AUTO LOAD 

6 

TIME STEP 

10.0 

CONTINUE 

0F xnput DECK’ 
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APPENDIX B 


Description of MARC Commands 
The MARC input deck consists of three blocks. They are as follows: 

Parameter Cards: 

TITLE: It gives the title given to the problem. This problem has the title as " /TEETH MODEL". 

SIZING: This specifics the size of the workspace buffer in number of words. This size should be 1.7 Mwords less 
than that specified in the JCL file. This is required for MARC program to run. Here, it is quite a large value. 

ELEMENT: This gives the element type used in the model. The one used here is number 7. 

PRINT, 5: It is a special printing option. ”5” means additional contact analysis information regarding nodes 
touching or separating from surfaces is given. 

END: It signals the end of Parameter Card block. 

Model Definition Cards: 

These cards contain the FE model data for the analysis. The data represents :- 

a) The FE mesh topology - element connectivity, nodal coordinates and sheet thickness. 

b) Material properties 

c) Loading and boundary conditions 

d) Nonlinear analysis control 

e) Output controls 

f) Contact analysis controls 
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FE mesh topology 

It is defined by the following cards:- 

CONNECnVITY : This defines connecdvuty for the elements in the model. A typical card is as illustrated : 
1 7 12345678 

First number denotes element number. 

“7* is the element type. It is HEX 8-noded clement. 

The last 8 numbers define the connectivity of the element in the counter clockwise direction. 

COORDINATES: This gives the coordinates for the nodes in the model. A typical card is: 

1 0.0 0.0 0.0 

"1” means 1st node number 

The other three numbers give the X, Y and Z coordinates of that particular node. 

The 2nd card is 

6 111 1 

"6" means max. number of coordinate directions to be read in per node. 

”111” means the total number nodes in the model 

"1" in the 20th column suppresses printout of the nodal coordinates in the output file. 

ISOTROPIC: This lets one define material properties, a yield criteria, etc. 

2nd card: 

no. of sets of isotropic material data to follow. 

3rd card: 
l.VON MISES 

”1” is material identification set no. 

"VON MISES” is the selected yield criteria 
4th card: 

3.00E+07,0.3 gives Young’s Modulus and poisson ratio 
5th card: 

1 TO # gives list of elements associated with this material 

SPRINGS: They are used to input any simple linear springs. 


90 



2nd card: 


103,1,51,1.50 

”103' gives the node to which 1st end of spring will be attached 
”1" gives the DOF at above node to which spring will be attached 
”51” gives node to which the other end of the spring will be attached 
”1” gives DOF at above node to which spring will be attached 
”50” gives the stiffness of the spring. 

Loading and Boundary Conditions 

FIXED DISP: It is used to prescribe displacement boundary conditions. 

2nd Card: No. of sets of boudary condition cards to be defined. 

3rd Card: 

"0.0, 0.0, 0.0" give the prescribed displacements for 1st, 2nd and 3rd DOFs 
i.e X, Y and Z directions. 

4th Card: 

1,2 

”1* means the X direction DOF 
”2” means the Y direction DOF 
5th Card: 

76,78,82 to 200 by 5.. 

This gives the list of nodes to which above boundary condition are applied. 

These above 3 Cards are repeated as required to define displacement at various nodes. 

POINT LOAD: 

2nd Card: no. of sets of point loads to be entered. It can be left blank. 

Here, it is left blank. 

3rd Card: 

0.0,10.0,0.0 

"0.0” gives nodal load associated with 1st DOF i.e X direction 
"10.0" gives nodal load associated with 2nd DOF i.e Y direction 
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"0.0" gives nodal load associated with 3rd DOF i.e Z direction 
4th Card: Gives list of nodes having the point load given above. 

"5141" is the node at which load is applied. 

Non-linear Analysis Control 

CONTROL: It lets one input parameters which control the convergence and accuracy of the non-linear analysis. 
2nd Card: 

204,0,1,0, 

”20” means max. no. of load steps 

"S" means max. no. of recycles during an increment 

"0" min. no. of recycles during an increment 

"1" this flags the convergence testing on displacements 

"0.0” flags for relative error testing 

"blank” default Full Newton Raphson iterative scheme is used 
3rd Card: 

".15" gives a relative error of 15% 

Output Controls 

POST: It creates a post processing tape for PATRAN. 

2nd Card: 

7 11 1 

"7* is no. of element variable to be written 

"1" is in 16-20 column for formatted POST TAPE 

"1” in 20-25 column is to write connectivity and coordinates on POST TAPE 

”1” in 45 column is to write post data every increment 

3rd Card: Gives various post codes 

11-16 give components of generalised stresses 

17 gives Equivalent Von Mises stress 

PRINT ELEM: Gives elements at which output is to be printed 
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2nd Card: 


1,2 

*1" is no. of sets 

"2" is increment between printout. Default is every increment printed. 

3rd Card: 

"STRAIN STRESS" are values to be printed 
4th Card: 

1 TO 44 is the list of elements to be printed 
5th Card: it lists integration points. 

PRINT NODE: gives information on nodal printout 
2nd Card: 

1,2 

"1" is no. of sets 

"2* is increment between printout. Default is every increment printed. 

3rd Card: 

”TOTA,LOAD,STRESS" are values to be printed 
4th Card: 

50 TO 85 is the list of nodes to be printed 

ERROR ESTIMATE: It gives error associated with FE discretization. Large values indicate stress gradients arc 
not accurately represented. 

2nd Card: 

1,1 

"1" in 1-5 column is for stress measure to be evaluated 
"1" in 6-10 column is for geometric measure to be evaluated 

SUMMARY: Gives the summary of output 

Contact analysis controls 
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CONTACT: Allows one to perform automated contact analysis without use of GAP elements for rigid to 
deformable contact as well as deformable to deformable contact. Here, it is deformable-deformable contact. 
2nd Card: 

2,60,60 

"2” tells two surfaces (bodies) will be defined 
"60* shows there are max. of 60 entities to be created for any body 
"60" shows max. no. of nodes that lie on the deformable surface 
3rd Card: 

,0.001, gives the distance below which node is considered touching a surface 
04th Card: 

1,0 

"1" 1st body 

"0" implies deformable body 
5th Card and 6th Card: 
blank 
7th Card: 

17 TO 32 gives list of elements for body one. 

8th TO 11th Card are repeated as above from 6th Card 

CONTACT TABLE: It is used for deactivating or activating bodies when the CONTACT option is used. 

2nd Card: 

"1" no. of sets of bodies to be input 
3rd Card: 

"1" gives the touching body number 
4th Card: 

"2" list of bodies for which the above body will detect contact 

OPTIMIZE: It allows a choice of bandwidth optimizers to be invoked . Is helpful in reducing the bandwidth 
1st Card: 

OPTIMIZER,,, 1 
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”4’ shows that Wavefront algorithm based on connectivity followed by Grooms algorithm is used. It is an 
effective technique for 3-D (complex) meshes. 

2nd Card: 

1 TO 44 gives list of elements 
3rd Card: 

”3* gives acceptable half-bandwidth at which it should exit down the iteration loop. 

Load Incrementation Cards: 

TIME STEP: Allows user to prescribe time step for static analysis 
2nd Card: 

*10” means time step of 10 sec 

AUTO LOAD: It describes number of equal load steps applied 
2nd Card: 

*2” shows 2 equal load increments 

CONTINUE: This terminates Load Incrementation or History Definition Cards. 
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APPENDIX C 


JCL AND USER SUBROUTINE 


The JCL file 


# USER-tobibel PW-tobibel 
I QSUB -r Xt7tr0 

# QSUB -IT 10800 
» QSUB -1M 30.2MW 

♦ QSUB -eo 

♦ QSUB -q systems 
set -x 
ja 

news 

dir-/hogs/tobibel/7teeth 

job-7tr0 
cd Sdir 

rm $dit/S job. 1st Sdir/S job. post deletes old output and post files Detore beginning new run 

Is -alp 

marc-/wrk/wmarc/marck5/marc defines the directory where MARC is accessed from 

S 

SAccesses MARC; defines input file, output file, and post file 
SIf RESTART option is used, then output restart file is written 
Sand also for subsequent restart the input restart file is mentioned 
S 

Snare i-S job. marc o-Sjob.lst patran-S job. post orestart-S job. restl usub-Sjob.f news-no 
Is -alp 
ja -st 


requester's name could be different from user 
CPU time requested 
memory size requested 

written only for priority systems queue 
gives news about MARC 

defines where the input file is present 
shows the job name 

changes to the directory where input file is present 


SUBROUTINE INLIST 


MARC USER SUBROUTINE TO DISABLE ECHO PRINTING OF THE INPUT DATA. 


C 

C 


WRITE (6,6000) 
RETURN 


6000 FORMAT ( 2 (/),' ENTRY TO * INLIST" USER SUBROUTINE’ ,6X,90( ’-•) , 

G 2(/),SX, ’ECHO OF THE INPUT DATA HAS BEEN DISABLED’, 

6 2 (/) , • EXIT FROM -INLIST” USER SUBROUTINE’ ,5X,90( •-•) ,/ ) 

END 
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MARC-CRAY K5-1-L 


NORTH AMERICAN MARC ANALYSIS RESEARCH CORPORATION 

EUROPE 


PACIFIC RIM 


IMPORTANT OUTPUT LISTING 



TRY TO "INLIST" USER SUBROUTINE ====== 

ECHO OF THE INPUT DATA HAS BEEN DISABLED 
IT FROM "INLIST" USER SUBROUTINE ====== 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

program sizing and options requested as follows 


element type requestedXXXXXXXXXXXXXXXXXXXXXXXXX 7 

number of elements in meshxxxxxxxxxxxxxxxxxxxxx 130A 
number of nodes in meshxxxxxxxxxxxxxxxxxxxxxxxx 1980 
max number of elements in any dist load listxxx 0 
maximum number of boundary conditionsxxxxxxxxxx 35 
load correction flagged or setxxxxxxxxxxxxxxxxx 
number of lists of distributed loadsXXXXXXXXXXX 3 

element out of core storage flagged, buffer size A096 

option for debug print outxxxxxxxxxxxxxxxxxxxxx 1 

stresses stored at all integration pointsxxxxxx 
tape no. for input of coordinates + connectivity 5 

no. of simple linear springs XXXXXXXXXXXXXXXXXX 12 

no, of different materials 2 max.no of slopes 5 

maximum elements variables par point on post tp 33 

number of points on shell section xxiixxxxxxxxxx 11 

option for terminal debugxxxxxxxxxxxxxxxxsxxxxx 
new style input format will be usedxxxxxxxxxxxx 
maximum number of set names isXXXXXXXXXXXXXXXXX 10 


number 

vector 


f processors used xxxxxxxxxxxxxxxxxxxxx 
i ^"^ukxxxxxxxkxxxxvXXXXXXXX 


c used 

xxxxxxxxX¥XXXXXXXXXXXXXXXXXX 


1 

1 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



key to stress, strain and displacement output 


element type 7 


8-node isoparametric brick 

stresses and strains in global directions 
l=xx 
2=yy 
3=zz 
4=xy 
5=yz 
6=xz 


displacements in global directions 
l=u global x direction 
2=v global y direction 
3=w global z direction 


KO 

'O workspace needed for input and stiffness assembly 191452 

i 

internal core allocation parameters 

degrees of freedom per node (ndeg) 3 

coords per node (nerd) 3 

strains per integration point (ngens) 6 

max. nodes per element (nnodmx) 8 

max. stress components per int. point (nstrmx) 6 

max. invariants per int. points (neqst) 1 


flag for element storage (ielsto) 1 

elems out of core, words per elem (nelsto) 1302 

elems per buffer (mxels) 3 

out-of-core space needed for element storage - 1699110 based on record size of 3906 

vectors in core, total space required 93060 


words per track on disk set to 4096 


internal element variables 
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internal element number 1 library code type 7 
number of nodes 3 8 

stresses stored per integration point = 6 

direct continuum components stored = 3 

shear continuum components stored = 3 

shell/beam flag = 0 

curvilinear coord, flag = 0 

int. points for elem. stiffness 8 

number of local inertia directions 3 

int. point for print if all points not flagged 9 

int. points for dist. surface loads (pressure) 6 

library code type = 7 

no local rotation flag = 1 

generalized displ. flag = 0 

large disp. row counts 6 6 6 11 11 11 


residual load correction is invoked 


$ XXX control cards for post-processor tape XXX 
$ specify post codes and labels for each element variable. 

$...! 1 ! 2 ! 3 ! 6 ! 5 f 6 • 7 ! 8 

post 


XXX note - format of post code cards has changed. 

in k6, enter code in first field and layer number in second field 

elem vars,post tape.prev tape, type , conn fl ,post tape, prev tape, repost , frequency, k2post 
6 16 17 1 1 19 20 0 3 0 

element variables appear on post-processor tape 16 in following order 
post variable 1 is post code 17 «= von mises stress 

post variable 2 is post code 131 = maximum principal stress 

post variable 3 is post code 133 ■ minimum principal stress 

6 XXX model definition cards Xxx 

$ specify components for each load and the nodes where it is applied. 

$...! . .. 1 ....!.... 2 ....!.... 3 .... I ....*....!.... 5 .... 7....!.... 8 

Spoint load 
$ 1 

$ 0.0 10.000 0.0 0.0 0.0 0 0 

$1826 

$ - specify displacements and the associated dofs and nodes for each set. 

2 . 3 ... 6 .. . 5 . 6 . .. 7. .8 


xxxmaximum record length on formatted post file= 


80 


approximate no. of records per increment on file 3 
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fixed disp 


fixed displacement » O.OOOE+OO O.OOOE+OO 
a list of degrees of freedom given below 
3 


a list of nodes given below 

905 1355 1990 1600 1609 1620 

fixed displacement * O.OOOE+OO O.OOOE+OO 


O.OOOE+OO 

1789 1800 

O.OOOE+OO 


a list of degrees of freedom given below 
12 3 

a list of nodes given below 
1976 

fixed displacement 3 O.OOOE+OO O.OOOE+OO O.OOOE+OO 
a list of degrees of freedom given below 
12 3 

a list of nodes given below 

15 955 520 700 729 790 889 900 


fixed boundary condition summary. 

total fixed degrees of freedom read so far 3 35 


b.c. 

node 

degree of 

magnitude 

b.c. 

node 

degree of 

magnitude 

number 


freedom 


number 


freedom 


1 

905 

3 

O.OOOE+OO 

2 

1355 

3 

O.OOOE+OO 

3 

1990 

3 

O.OOOE+OO 

9 

1600 

3 

O.OOOE+OO 

5 

1609 

3 

ii OOOE+OO 

6 

1620 

3 

0.000E+00 

7 

1789 

3 

0.000E+00 

8 

1800 

3 

0.000E+00 

9 

1976 

1 

O.OOOE+OO 

10 

1976 

2 

O.OOOE+OO 

11 

1976 

3 

0.000E+00 

12 

15 

1 

O.OOOE+OO 

13 

15 

2 

0.000E+00 

19 

15 

3 

O.OOOE+OO 

15 

955 

1 

0.000E+00 

16 

955 

2 

0.000E+00 

17 

955 

3 

O.OOOE+OO 

18 

520 

1 

O.OOOE+OO 

19 

520 

2 

0.000E+00 

20 

520 

3 

O.OOOE+OO 

21 

700 

1 

O.OOOE+OO 

22 

700 

2 

O.OOOE+OO 

23 

700 

3 

0.000E+00 

29 

729 

1 

O.OOOE+OO 

25 

729 

2 

O.OOOE+OO 

26 

729 

3 

0.000E+00 

27 

790 

1 

O.OOOE+OO 

28 

790 

2 

O.OOOE+OO 

29 

790 

3 

0.000E+00 

30 

889 

1 

O.OOOE+OO 

31 

889 

2 

O.OOOE+OO 

32 

889 

3 

O.OOOE+OO 

33 

900 

1 

0.000E+00 

39 

900 

2 

O.OOOE+OO 

35 

900 

3 

0 . 900E+00 





< --- 

specify 
. . .1 

connectivity 
I 2 1 . 

for each element. 

. . .5. . . . ? « . I 

_ . 1 _ 


.1 .... 7 1 

A 
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connectivity 


meshrl,iprnt 
5 1 

elem no., type, nodes 

* specify coordinates for each node. 

»...! 1....I 2 I....3 ! 4. ... I ... .5. ... ! 6. ... I 7 ! B 

coordinates 


ncrdl ,meshrl,iprnt 
3 5 1 

node coordinates 

isotropic 


isotropic material material id ■ 1 

von mises yield criteria 
isotropic hardening rule 

a nu rho alpha yield yield2 

0 . 300E+0B 0 . 300E+00 O.OOOE+OO O.OOOE+OO 0.100E+21 0.100E+21 

from element 1 to element 1224 by 1 


isotropic material material id ■ 2 

von mises yield criteria 
isotropic hardening rule 

e nu rho alpha yield yield2 

0.300E+11 0 . 300E+00 O.OOOE+OO O.OOOE+OO 0.100E+21 0.100E+21 

from element 1225 to element 1304 by 1 


springs 


spring 

node 

degree 

node 

degree 

spring damping 

number 


freedom 


freedom 

force 

force 

1 

11 

1 

995 

1 

1 . 000E+02 

0.000E+00 

2 

11 

2 

995 

2 

1 . OOOE+02 

O.OOOE+OO 

3 

11 

3 

995 

3 

1 . 000E+02 

O.OOOE+OO 

4 

96 

1 

905 

1 

1. 0O0E+O2 

O.OOOE+OO 

5 

96 

2 

905 

2 

1. 000E+02 

0.000E+00 

6 

96 

3 

905 

3 

1. 000E+02 

0.000E+00 

7 

451 

1 

1400 

1 

1. 000E+02 

O.OOOE+OO 

8 

451 

2 

1400 

2 

1. 000E+02 

0.000E+00 

9 

451 

3 

1400 

3 

1. 000E+02 

O.OOOE+OO 

10 

496 

1 

1355 

1 

1. 000E+02 

O.OOOE+OO 

11 

496 

2 

1355 

2 

1. OOOE+02 

O.OOOE+OO 

12 

496 

3 

1355 

3 

1. 000E+02 

O.OOOE+OO 
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contact 


number of bodiaa ■ 2 
max number of entities per body ■ 250 
bound on number of boundary nodes * 250 
friction type(l~m , 2-coulomb) = 0 
distrib-0 or nodal-1 coul. frict “ 0 


relative velocity below which a 

node is considered sticking 3 O.OOOOOE+OO 

distance below which a node is 

considered touching a surface “ O.OOOOOE+OO 

nodal reaction above which a node 

separates from a body a O.OOOOOE+OO 


body number = 1 

number of sets of data = 0 


body positioning data 
1st coordinate of center of rotation 
2nd coordinata of canter of rotation 
3rd coordinate of center of rotation 
1st component of velocity 
2nd component of velocity 
3rd component of velocity 


0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 


body positioning data continue 
angular velocity 0.00000 

total angle rotated around axis 0.00000 

1st component of directional cosine 0.00000 

2nd component of directional cosine 0.00000 

3rd component of directional cosine 0.00000 

friction coefficient 0.00000 

from element 616 to element 936 by 6 


body number ■ 2 

number of sets of data m 0 


body positioning data 

1st coordinate of center of rotation 0.00000 

2nd coordinate of center of rotation 0.00000 

3rd coordinate of center of rotation 0.00000 

1st component of velocity 0.00000 

2nd component of velocity 0.00000 

3rd component of velocity 0.00000 

body positioning data continue 
angular velocity 0.00000 

total angle rotated around axis 0.00000 

1st component of directional cosine 0.00000 

2nd component of directional cosine 0.00000 

3rd component of directional cosine 0.00000 

friction coefficient 0.00000 
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1 to •lament 321 by 


4 


from element 
and 

a list of elements given below 
325 341 351 373 389 

421 437 597 


405 


contact table 


2 

flexible body number 1 and flexible body number 2 will be detected by each other 
optimized , , , ,1 


cuthill-mckee algorithm 
control 


max. max. min. 

incs recycles recycles 

240 5 0 

maximum allowed relative error in residual forces 0 . 15000E+00 


full newton- raphson technique chosen 
restart 

iwhich<incbrs,incrs,out no. ,in no. »itri,ireprt,ilasti,incsur .ilasts 
1408 900000 

print elem 


values will be printed at integration points 
element quantities printed every 4 increments 
strain.stress 

from element 616 to element 936 by 4 

and 

from element 1 to element 321 by 4 

and 

a list of elements given below 
325 341 351 373 

421 437 597 


389 


405 
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print node 


number of sets used for selective print of nodal quantities is 1 
nodal quantities printed every 4 increments 
reac , to ta , stress 

from node 1 to node 491 by 10 

and 

from node 6 to node 496 by 10 

and 

from node 701 to node 781 by 20 

and 

from node 801 to node 881 by 20 

and 

from node 905 to node 1400 by 5 

error estimate 


stress discontinuity will be evaluated 
geometric distortion will be evaluated 

end option 


maximum connectivity is 18 at node 745 


workspace needed for optimizing ■ 248979 


maximum connectivity is 24 at noda 1980 


maximum half-bandwidth is 222 between nodes 1500 and 1721 
number of profile entries including fill-in is 224672 

number of profile entries excluding fill-in is 21392 

XXXKXXXXXXX 

distance below which a node is considered 

touching a surface is 2.51180E-03 

xxxxxxxxxxx 


total workspace needed with in-core matrix storage ■ 2246772 



106 


restart information 


...file 8-- maximum record length 3 83160 

approximate no. of words on file per increment 3 


sliding velocity below which sticking is considered 0.10000E+05 


warning! input data was not enough for this value to 
results may be wrong if friction calculations are 


be calculated, 
being made. 


load increments associated with each degree of freedom 
summed over the whole model 


distributed loads 
0 . 000E+00 0 . 000E+00 Q.OOOE+OG 


point loads 

0 . OOOE+OO 0 . 000E+00 0. OOOE+OO 


increment zero is a null step 


worst original aspect ratio is 10.003 at element 32$ 
worst original warpage ratio is 2.050 at element 289 

worst current aspect ratio is 10.003 at element 329 
worst current warpage ratio is 2.050 at element 289 

largest change in aspect ratio is 1.000 at element 20 
largest change in warpage ratio is 1.000 at element 13 


generalized stresses 


1 0 . OOOOOE+OO 

6 0. OOOOOE+OO 

11 0. 00000E+00 


0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


0. OOOOOE+OO 
0. OOOOOE+OO 
O.OOOOOE+OO 


0. OOOOOE+OO 
O.OOOOOE+OO 
0.00000E+00 


O.OOOOOE+OO 

O.OOOOOE+OO 

O.OOOOOE+OO 


O.OOOOOE+OO 

0.00000E+00 

0.00000E+00 


1867308 


T«» 
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forcts in linear springs 
spring no. nodal dofl nodo2 dof2 force 


1 

11 

1 

995 

1 

0 . 0000E+00 

2 

11 

2 

995 

2 

0 . OOOOE+OO 

3 

11 

3 

995 

3 

0 . 0D00E+00 

9 

96 

1 

905 

1 

0 . 0000E+00 

5 

96 

2 

905 

2 

0 . OOOOE+OO 

6 

96 

3 

905 

3 

0. OOOOE+OO 

7 

951 

1 

1900 

1 

0. OOOOE+OO 

8 

951 

2 

1900 

2 

0. OOOOE+OO 

9 

951 

3 

1900 

3 

0. OOOOE+OO 

10 

996 

1 

1355 

1 

0. OOOOE+OO 

11 

996 

2 

1355 

2 

0. OOOOE+OO 

12 

996 

3 

1355 

3 

0. OOOOE+OO 


end of increment 0 

restart date at increment 0. 0 on tape 8 

formatted post data at increment 0. 0 on tape 19 

time * 15.67 

point load 


read from unit 5 
0 . 00DE+00 0.600E+01 O.OOOE+OO 
a list of nodes given below 
1732 

auto load 


iotnum,incasm 
1 0 

time step 


time increment “ 10.00000 

continue 


equal load incs specified for 1 increments 
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start of incromont 1 


load increments associated with each degrea of freedom 
summed over the whole model 

distributed loads 
O.OOOE+OO 0 . 000E+00 0.000E+00 


point loads 

O.OOOE+OO 6 . 000E+00 O.OOOE+OO 

start of assembly 
time * 16.04 


start of matrix solution 
time ■ 21.51 


singularity ratio 6.3282E-04 


end of matrix solution 
time = 33.11 

maximum residual force at node 1905 degree of freedom 2 is equal to 
maximum reaction force at node 455 degree of freedom 2 is equal to 
convergence ratio 


0 . 112E-04 
0 . 263E+01 
0.426E-05 


separation force required is 0.11210E-04 


MARC-CRAY K5-1, 06/13/93, output for increment 


1. with 2 prop. ids 


total transient time 


1 . 00000E+01 



109 


worst current aspect ratio is 10.004 at element 324 
worst current warpage ratio is 2.050 at element 289 

largest change in aspect ratio is 1.000 at element 936 
largest change in warpage ratio is 1.000 at element 940 


largest normalized stress jump is 0.79283E+12 at node 1442 component 


4 mean value is 


largest stress jump is 0.25036E+04 at node 479 component 1 mean value is - 

end of increment 1 
time = 41.53 

point load 


read from unit 5 
0 . 000E+00 0 . 200E+01 0.000E+00 
a list of nodes given below 
1732 

auto load 


iotnumjincasm 
3 0 

time step 

time increment “ 10.00000 

continue 


equal load incs specified for 3 increments 


start of i 


node 1260 is touching body 2 patch 149 
the patch it touches is 149 
the internal node numbers for it are 
the normal vectors are -0.73027 

node 1210 is touching body 2 patch 126 
the patch it touches is 126 
the internal node numbers for it are 
the normal vectors are -0.70847 


n c r e 

men 

t 

2 

391 

396 

346 

341 

0.67742 


-0.00839 

336 

341 

291 

286 

0.69556 


-0.11947 


0.48236E-10 
0 . 10067E+04 
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3 

11 

4 

96 

5 

96 

6 

96 

7 

451 

8 

451 

9 

451 

10 

496 

11 

496 

12 

496 


3 

995 

3 

-0.1869 

1 

905 

1 

0.2886 

2 

905 

2 

1.800 

3 

905 

3 

-6 . 9383E-03 

1 

1400 

1 

-2.7481E-02 

2 

1400 

2 

2.121 

3 

1400 

3 

-6.1306E-02 

1 

1355 

1 

-0.1154 

2 

1355 

2 

2.109 

3 

1355 

3 

-0.1131 

i d 

o f 

i n 

c r e m e n t 


restart data at increment 8. 0 on tape 8 

time ■ 352.85 


**# end of input deck - job ends 


marc exit number 3004 
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